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ABSTRACT 


This  thesis  explores  a  MATLAB  implementation  of  a  Fourier  transform  approach 
to  model  and  predict  transient  optical  wave  propagation  through  free-space.  A  three-step 
approach  is  adopted  in  this  study.  First,  the  mathematical  development  establishes  the 
importance  of  the  total  impulse  response  as  the  Green's  function,  meeting  the  boundary 
conditions  and  solving  the  wave  equation.  Second,  a  MATLAB  program  is  developed  to 
simulate  the  mathematical  model  by  computing  and  displaying  the  graphical 
representation  of  an  optical  wave's  spatial  distribution  on  a  plane  at  a  given  distance  from 
a  spatially  filtered  source.  Third,  a  circular  excitation  function  is  used  to  verify  the 
program  and  then  the  results  of  another  three  excitations,  namely  the  square,  circularly 
truncated  Gaussian  and  circularly  truncated  Bessel  functions  are  similarly  generated.  The 
effort  of  this  thesis  provides  an  inexpensive  means  to  analyze  a  transient  optical  wave 
propagation  of  a  spatially  filtered  optical  source. 
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1.  INTRODUCTION 


Advances  in  laser  technology  have  made  coherent  optical  sources  readily 
available.  With  applications  such  as  image  processing,  image  pattern  recognition, 
spectrum  analysis,  synthetic-aperture  radar  data  processing  and  biomedical  applications, 
laser  sources  may  be  broadly  classified  into  continuous  and  pulsed  lasers.  Any  laser 
which  operates  for  a  second  or  more  at  a  time  is  called  "continuous  wave."  There  are  also 
many  other  types  of  lasers  that  operate  only  in  the  pulsed  mode.  For  example,  in  solid- 
state  lasers,  the  key  problem  is  heat  dissipation.  It  takes  time  for  excess  pump  energy 
delivered  to  the  laser  rod  to  make  its  way  out  as  heat  and  continuous  wave  operation  can 
cause  heat  to  build  up  to  laser  damaging  levels. 

The  pulsed  mode  laser  also  finds  many  other  applications  which  exploit  its  short 
pulse  duration.  The  short  length  of  the  pulse  makes  it  an  ideal  candidate  for  three- 
dimensional  imaging,  either  to  acquire  depth  resolution  through  range  gating  or  to 
discriminate  against  scattering.  Very  high  peak  intensities  can  be  reached  at  moderate 
pulse  energies  with  ultrashort  pulses.  Finally,  the  ability  to  make  nondispersive  "solitons" 
led  to  a  new  pulsed  code  communication  system  with  optical  fibers.  All  these  desirable 
qualities  of  a  laser  source  are  made  possible  solely  because  of  the  short  pulse  duration 
and  this  continuing  exploitation  has  led  scientists  to  discover  new  techniques  to  produce 
ultrashort  pulses  in  the  femtosecond  regime.  This  thesis  tries  to  find  a  method  to  model 
and  predict  the  behavior  of  laser  pulse  propagation  using  computer  simulation. 
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A. 


PROPAGATION  OF  PULSED  FIELDS 


Laser  sources  exhibit  a  spatial  amplitude  distribution,  which  is  typically  Gaussian. 
It  is  possible  to  spatially  filter  such  a  beam  to  produce  an  alternative  shaped  beam.  Such  a 
variation  may  exhibit  a  circular  or  square  uniform  cross-section  and  either  of  these  could 
have  an  arbitrary  spatial  weighting  distribution.  The  utility  of  such  filtering  is  unknown 
unless  the  diffracted  field  distribution  can  be  predicted  at  any  given  distance. 

The  theory  of  linear  systems  can  be  applied  for  our  purpose  of  predicting  this 
diffracted  field  distribution.  By  taking  the  multi-dimensional  Fourier  transform  of  the 
complex  field  distribution  across  any  plane,  the  spatial  Fourier  components  can  be 
identified  as  plane  waves  travelling  in  different  directions.  Accounting  for  phase  shift 
during  travel  and  applying  the  superposition  theorem,  the  field  amplitude  at  any  other 
point  will  be  the  sum  of  each  of  these  contributing  waves  directions.  Thus,  the 
propagation  phenomenon  of  the  optical  wave  may  be  regarded  as  a  linear  space-invariant 
system  characterized  by  a  specific  transfer  function. 

B.  PROBLEM  DESCRIPTION 

In  this  thesis,  we  want  to  consider  the  prediction  of  transient  optical  waves  after 
free-space  propagation  from  one  plane,  where  the  wave  is  known,  to  a  parallel  plane  that 
is  located  a  distance  z  away.  We  shall  denote  these  two  planes  as  the  input  plane  and 
output  plane,  respectively.  The  assumed  geometry  is  shown  in  Figure  1. 

The  wave  is  assumed  known  in  the  z  =  0  plane  and  is  given  by  uo(x,y,0,t).  For  our 
model,  the  input  excitation  must  be  separable  in  space  and  time,  i.e., 

Ug(x,y,0,t)  =  s(x,y,0)T(t)  (1) 
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where  s(x,y,0)  is  the  spatial  portion  of  the  excitation  and  T(t)  is  the  temporal  portion  of 
the  excitation.  The  propagation  medium  is  assumed  linear  and  homogeneous;  in  this 
thesis,  we  assume  free  space. 


output  field,  u(x,y,z,t) 


input  field,  uo{x,y,0,t) 


Figure  1.  Assumed  geometry 


We  use  scalar  wave  theory  [Ref  1]  to  represent  the  optical  wave.  Our  aim  is  to 
predict  u(x,y,z,t)  on  the  oufr)ut  plane,  given  uo(x,y,0,t)  on  the  input  plane  and  the  distance 
z  unit  away  from  the  source  plane.  The  constraints  are  that  the  wave  must  solve  the  scalar 
wave  equation 


vv2  .  .X  1  ^  u(x,y,z,t)  ^ 

V  u{x,y,z,t)-— - - =  0 

c  d  t 


(2) 


and,  since  we  are  considering  propagation  in  free-space  (i.e.,  no  boundaries  are  present 


other  than  at  the  input  plane),  the  wave  goes  to  zero  as  the  distance  r  =  -yjx^  +y^  +z'^ 
goes  to  infinity  in  the  half-space  above  of  the  input  plane. 
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In  addition,  we  have  also  made  the  following  assumptions  that  helped  to  simplily 
our  study  and  simulation  of  the  field  distribution  in  the  output  plane: 

1.  We  have  fixed  the  size  of  the  input  and  output  planes,  so  that  we  may 
concentrate  our  observation  and  analysis  on  the  center  area  of  the  wave 
distribution  on  these  planes. 

2.  We  accounted  for  the  effects  of  diffraction  by  using  suitable  Green's 
fimction  for  our  model. 

3.  We  fixed  the  distance  z,  between  the  source  and  image  planes  so  that  we 
may  plot  u(x,y,z,t)  in  three-dimensional  graphical  representations. 

4.  We  considered  variable  aperture  sizes  to  suit  different  input  excitations. 
The  approach  adopted  to  solve  our  thesis  propagation  question  may  be 

summarized  into  a  flow  chart  as  shown  in  Figure  2.  A  mathematical  expression  based  on 
linear  system  theory  and  Fourier  transform  is  derived  for  the  predicted  field,  u(x,y,z,t). 
This  expression  is  then  developed  into  a  MATLAB  program,  which,  given  a  knovra 
excitation  at  the  input  plane,  predicts  (or  simulate)  the  expected  field  distribution  at  the 
output  plane. 

Four  input  excitation  functions  were  used: 

1 .  Circular  field  distribution, 

2.  Square  field  distribution, 

3.  Circularly  truncated  Gaussian  field  distribution  and 

4.  Circularly  truncated  Bessel  field  distribution. 
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We  shall  step  through  the  process  of  generating  the  predicted  field  distribution  for 
a  circular  input  excitation  in  order  to  verify  the  accuracy  of  our  MATLAB  programs  and 
then  results  of  the  other  three  excitations  will  be  generated. 


Figure  2.  Research  approach  for  thesis. 

Now  that  we  have  defined  the  scope  of  our  thesis  research,  next  we  shall  discuss 
the  theories  involved. 


5 


6 


II.  THEORY 


Two  main  theories  were  involved  in  this  thesis  research:  linear  systems  theory  and 
Fourier  transform  theory.  Section  A  develops  the  concept  of  how  linear  system  theory 
may  be  used  to  characterize  wave  propagation  model  in  terms  of  a  transfer  function  (also 
known  as  the  spatial  impulse  response  or  Green’s  function).  Section  B  shows  how  the 
field  distribution  at  the  output  plane  may  be  found  by  solving  the  wave  equation  using  a 
set  of  defined  propagation  and  boimdary  conditions  specified  in  our  problem  description 
in  Chapter  I.  Section  C  shows  that  the  temporal  spatial  impulse  response  may  be  derived 
from  the  expression  of  the  computed  field  distribution  at  the  output  plane.  Section  D 
demonstrates  how  the  temporal  spatial  impulse  response  may  be  expressed  in  a  suitable 
form  for  computer  simulation  by  taking  its  spatial  Fourier  transform.  Finally,  Section  E 
provides  an  overview  of  the  software  that  was  used  in  our  simulation  program. 

A.  PROPAGATION  MODEL  AS  A  LINEAR  SYSTEM 

Many  physical  phenomena  are  found  to  share  the  basic  property  that  their 
response  to  several  stimuli  acting  simultaneously  is  identically  equal  to  the  sum  of  the 
responses  that  each  stimulus  would  produce  individually.  Such  phenomena  are  called 
linear  and  the  property  they  share  is  called  linearity.  Optical  propagation  in  linear 
homogeneous  media  is  such  a  phenomenon.  The  wave  equation  (Equation  2)  leads  us  to 
regard  optical  propagation  as  a  linear  mapping  of  the  input  light  distribution  into  the 
output  light  distribution.  Therefore  we  may  consider  the  mapping  of  wave  distribution, 
uo(x,y,0,t)  to  u(x,y,z,t)  on  a  plane  located  z  unit  distance  away  as  linear  and  apply  all  of 
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the  properties  of  linear  system  in  simplifying  the  mathematics  that  describe  this 
operation. 

In  linear  system  theory,  we  characterize  a  mapping  operation  by  its  impulse 
response.  As  shown  in  Figure  3a,  the  impulse  response,  h(x,y,z,t),  is  the  response  of  the 
operation  to  an  input  of  S(x,y,z,t)=S(x,y,z)S(t).  In  propagation  terms,  the  impulse  response 
is  called  the  Green's  function,  which  is  the  solution  of  the  wave  equation  and  its  boundary 
conditions  to  an  impulse  excitation. 


5{x,y,z)5{t) 


Propagation  & 
boundary  conditions 


(a) 


s(x,y,0)5(t 


Propagation  & 
boundary  conditions 


(b) 


h(x,y,z,t) 

- > 


dh(x,y,z,t) 

p(x,y,  2,  t)  =  -s( X,  y,0)  *  *  - - - - 

X  y  a  z 

— > 

\dh(x,y,z,t) 

p(x,y,z.t)  =  -s(f.f^,0)3\ - - - 


Figure  3.  (a)  Impulse  response  (or  Green's  function)  and  (b)  temporal  impulse  response. 

Also,  as  we  shall  see  in  section  B,  we  may  predict,  for  a  spatially  invariant 
system,  the  response  to  a  source  with  an  arbitrary  excitation  and  impulse  temporal 
excitation  in  terms  of  the  impulse  response.  As  shown  in  Figure  3b,  if  we  represent  the 
input  excitation  on  the  source  plane  as  s(x,y,0)5(t),  the  field  distribution  at  the  output 
p\dxit,p(x,y,z,t),  will  be  given  as 


p( X,  y,  z,  t)  =  -s( X,  y.O)  *  * 

jt  y 


dh(x,y,z,t) 

dz 


(3) 
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where  *  notation  indicates  convolution  over  the  variable  noted  [Ref.  2].  We  shall  call  the 
output  field  distribution,  p(x,y,z,t),  the  temporal  spatial  impulse  response  (i.e.,  it  is  the 
response  of  the  system  to  an  arbitrary  spatial  excitation  with  an  impulsive  temporal 
excitation). 

As  we  know  from  convolution  theory,  the  spatial  convolution  of  Equation  3  may 
be  converted  into  multiplication  in  the  spatial  frequency  domain  by  taking  its  spatial 
Fourier  transform  [Ref  5]; 


P(fx’fy’^’0  =  -s(fxJy>0)^< 


^h(f^,f,z,t) 


dz 


(4) 


Also  as  we  shall  see  later  that  for  our  computer  simulation  purposes.  Equation  4  is  a  more 
suitable  form  for  quick  computation  than  Equation  3. 

In  a  more  general  form,  the  output  field  distribution,  ^x,y,z,t),  to  an  excitation 
with  an  arbitrary  spatial  and  temporal  dependence  can  be  expressed  in  terms  of  the 
temporal  spatial  impulse  response  as  [Ref.  2] 


(l)(x,y,z,t)  =  T(t)  ^p(x,y,z,t).  (5) 

B.  SOLUTION  TO  WAVE  EQUATION 

To  derive  the  impulse  response,  we  first  need  to  find  a  solution  to  the  wave 
equation  meeting  the  set  of  propagation  and  boundary  conditions  defined  by  our 
propagation  model  in  Chapter  I.  From  [Ref  2],  the  solution  to  the  wave  equation. 
Equation  2,  is  given  by  the  radiation  integral.  Assuming  a  planar  input  aperture,  the  field 
u(x,y,z,t)  is  given  by 
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^Ug(x,y,0,t)  dh(x,y,z,t) 

u(  X,  y,  z,  t)  = - - - ***h(x,y,z,t)-Ug(x,  y,0,  t)*** - ^ — 

C/n  X  y  t  X  y  I 


(6) 


where  the  quantity  uo(x,y,0,t)  is  the  scalar  wave  distribution  at  the  source  plane,  h(x,y,z,t) 
is  the  Green's  function  that  both  solves  the  wave  equation  and  meets  the  boundary 
conditions  and  the  derivative  with  respect  to  n  represents  the  normal  derivative.  For  input 
and  output  planes  that  are  normal  to  the  z-axis,  the  normal  derivative  will  become  the 
derivative  with  respective  to  z.  Hence  Equation  6  may  be  rewritten  as 


dug(x,y,0,t)  dh(x,y,z,t) 

u(  X,  y,  z,  t)  = - - - ***h(x,y,z,t)-Ug(x,  y,0,  t)*** - - - .  (7) 

C/Z  ^  y  I  X  y  I  g  z 


In  this  thesis,  the  value  of  the  field  on  the  planar  source  plane,  uo(x,y,0,t),  is 
known.  Hence,  it  is  desirable  to  eliminate  the  normal  derivative  of  Equation  7  (i.e.,  the 
first  term  on  the  right  side  of  the  equation)  and  to  use  the  second  known  term  for  the 
solution.  This  can  be  done  by  using  a  Green's  function  given  by  [Ref  2]  which  has  also 
considered  the  effects  of  diffraction; 


10 


(10) 


dh  2zs(t--^'j  2zS'lt-^j 

dn  dz  cR^ 

where  R  =  yjr^  +z^  =  +y^  +z^  and  6'  indicates  the  time  derivative  of  the  Dirac 

delta  function. 

By  eliminating  the  known  first  term  on  the  right  of  Equation  7  and  substituting  the 
Green's  function  of  Equation  8,  the  field  can  then  be  written  as 

dh(x,y,z,t) 

u(  X,  y,  z,  t)  =  -Uq  ( X,  y,0,  t)*** - - - 

xy  t  dz 


=  Uo(x,y,0,t)***- 

X  y  t 


2zS 


,  +u,(x,y,0,t)***-  , 

K  X  y  I  CK 

This  equation  represents  the  expression  for  the  field  distribution  at  the  output  plane. 
C.  COMPUTATION  OF  THE  TEMPORAL  SPATIAL  RESPONSE 


i‘-yc) 


(11) 


To  simplify  Equation  1 1  further  so  that  it  is  easier  for  computer  simulation,  first 
we  take  its  two-dimensional  Fourier  transform  to  convert  convolution  in  the  space 
domain  into  multiplication  in  the  spatial  frequency  domain.  Then,  by  substituting 
Equation  1  into  the  expression,  we  have 


u(f^,fy,z,t)  =  3lu(x,y,z,t)}  =  2^Uo(x,y,0,t)*** 


dh(x,y,z,t) 


=  T(t)^ 


dz 

\dh(x,y,z,t) 


dz 


By  comparing  Equation  5  with  EquationlS,  we  observe  that 


(13) 
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(14) 


P(fxJy>^>t)  =  -s(f,,fy  ,0;3| 


\^h(x,y,z,t) 


dz 


The  inverse  spatial  transform  of  produces  the  spatial  impulse  response, 

p(x,y,z,t),  which  is  our  required  field  distribution  at  z  when  the  excitation  at  the  input 
plane  is  a  temporal  pulse  excitation. 

D.  TEMPORAL  SPATIAL  REPONSE  FOR  COMPUTER  SIMULATION 

We  now  want  to  find  an  expression  for  the  spatial  impulse  response,  p(x,y,z,t). 
Substituting  ( x,  y,0,  t)  =  s(x,y,0)d(t)  into  Equation  1 1 ,  we  have 


p(  X,  y,  z,  t)  =  s(  X,  y,0)S  (t)*** 


2z5[t-%) 


X  y  t 


+  s(x,y,0)d  (t)*** 

I  y  I 


2z5 


cR^ 


(15) 


Since  f  *  g'  =  (f  *  g)'  =  f '  *  g,  we  can  interchange  the  order  of  the  derivative  in  the  second 
term  of  Equation  15  and,  by  expanding  the  convolution  term  with  6(t),  get 
p(x,y,z,t) 


2z5 


('-%) 


=  s(x,y,0)6  (t)*** - - +  s(x,y,0)5  (t)*** 

X  y  !  K  X  y  t 


2z5 


{-%) 


cR^ 


(16) 


2zs[t-y^ 


=  s(x,y,0)** - '  ,  '  +S  (t)* 

X  y  K 


2zS{l-%) 


s(x,y.O)**-  , 

X  y  cR 


The  spatial  convolutions  over  x  and  y  in  last  line  of  Equation  16  are  more  easily 
performed  in  the  transform  domain.  Taking  the  two-dimensional  spatial  Fourier 
transform  of  Equation  16  gives 


12 


P(fx>fy>^>0  = 

's(f„fy,0)2zJo(pylch^  -z^)  (s(f„fy,0)2zJo(p4c^?'^j]  (l?) 

:je  - 

cV  ,1^  J 

where  p  is  the  two-dimensional  spatial  transform  of  p,  fx  and  fy  are  the  spatial 
frequencies,  p  is  the  radial  spatial  frequency  (p  =  ^fx^  +  //  )  and  the  transform  pan- 
given  below  has  been  used 


Y”  I  "  (ct)"-^  ^  ■ 


Recognizing  that  time  convolution  with  the  time  derivative  of  the  Dirac  impulse  is  the 
same  as  taking  the  derivative  of  the  function  in  the  time  domain,  i.e.. 


S  (t)*f(t)  =  f'(t), 


we  have 


p(x.y,z,t)  = 

+ Yt  ^  '^0 (pVcV-z^)/f(r  -  z/c)| 
By  factoring  the  common  term  's(f^,fy,,0)  from  Equation  20,  we  have 

p( X, y, z. t)  =  3'^ I? r/,, 4 ,0)  Jo )if(/  - zjc) 


which  may  be  simplified  further  as 
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(22) 


p(  X,  y,  z,  t)  =  3'^  \s(f^,  fy  ,0)c  ( /, ,  fyZ,  /;} 
where  we  denote  the  following  as  the  propagation  spatial  filter  [Ref.  2] 


For  our  simulation  model,  we  have  evaluated  c(f^,fy,z,t)  for  three  different  time 


regions, 


where  we  have  made  the  assumptions  for  the  second  line  of  Equation  24  that  both 
H{t  -  z/c)  and  S{t  -  z/c)  are  equal  to  one. 

Equations  22  and  24  are  the  only  two  equations  required  for  our  simulation 
program.  Recall  that  's(f^,fy,0)  represents  our  spatial  pulse  excitation  at  the  input  plane, 
p(f^,fy,z,t)  represents  the  spatial  field  distribution  at  the  output  plane  and 
c(f^,fy,z,t)  represents  the  linear  system  transfer  function  that  maps  's(f^,fy,0)  onto 
p(f^,fy,z,t).  In  optics  term,  c( f^,fy,z,t)  is  also  known  as  the  propagation  spatial 
filter.  It  is  the  "filtering  function"  that  meets  the  set  of  defined  boundary  conditions  in 
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Chapter  I.  It  also  characterizes  the  effects  of  diffraction  that  modify  the  input  excitation 
as  it  propagates  through  the  free  space  between  the  input  and  output  planes. 

E.  MATLAB  OVERVIEW 

MATLAB  is  an  acronym  for  MATrix  LABoratory.  It  is  a  high-performance, 
interactive,  scientific  and  engineering  software  package.  As  its  name  suggests,  its  basic 
data  element  is  a  matrix.  A  major  advantage  of  MATLAB  is  that  traditional  programming 
is  not  needed  since  problems  and  solutions  are  expressed  just  as  they  would  be  written 
mathematically.  Another  distinct  advantage  is  MATLAB's  expansion  capability  with 
preprogrammed  functions,  such  as  the  calculation  of  two-dimensional  FFTs  and  the 
calculation  of  Bessel  functions.  [Ref  7] 

There  are  two  types  of  macro-like  files  called  m-files  (called  m-files  for  the  ".m" 
suffix);  one  is  known  as  the  script  m-file  and  the  other  is  the  junction  m-file.  A  script  m- 
file  is  used  to  automate  long  sequences  of  commands  including  functions.  Arguments  are 
not  passed  into  script  files.  A  function  m-file,  however,  may  have  arguments  passed  into 
them.  Another  difference  between  the  two  file  types  is  that  the  first  line  of  a  function  m- 
file  begins  with  the  word  "function"  and  all  variables  used  in  the  function  are  local. 
Examples  of  script  m-files  in  this  thesis  are  IOPTFIL.m,  lOPTPROP.m, 
PLOTFILTER.m,  PLOTFIELD.m,  ANIMATEl.m,  ANIMATE2.m  and  ANIMATES .m 
(Appendixes  B,  C,  D  and  E).  Examples  of  function  m-files  include  the  input  excitation 
functions:  CRCLE.m,  SQUARE.m,  CRCGAUS.m  and  CRCBESS.m  (Appendix  A),  the 
three-dimensional  graphing  function  mesh  and  the  two  fft  functions  that  realize  the 
Fourier  transform  required  for  our  programs. 
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The  two  fft  functions  employed  for  Fourier  transform  are  fft2  and  fftshift;  fft2 
carries  out  a  two-dimensional  Fourier  transform  while  fftshift  carries  out  a  center-to- 
comer  geometrical  shift  on  an  input  function.  Both  of  these  functions  must  be  used 
together  and  they  perform  the  fast  Fourier  transform  on  an  input  function.  Since  these  are 
frequently  used  functions  throughout  our  program,  it  is  worth  the  attention  here  to 
elaborate  on  its  proper  usage,  especially  on  their  order  of  application  to  ensure  the  correct 
phase  result  is  obtained  for  the  resultant  function. 

The  correct  way  to  do  a  fast  Fourier  transform  in  MATLAB  is  to  do  it  in  three 
separate  operations.  First  a  fftshift  must  be  applied  to  the  input  function  which  we  shall 
denote  as  shft-input  and  its  result  as  input.  (The  "shft"  prefix  here  is  to  remind  us  that  a 
fftshift  operation  must  be  applied  first  prior  to  a  fft2  function.)  Then  a  fft2  function  is 
applied  and  the  result  is  denoted  as  F-input.  Finally  another  fftshift  function  is  applied 
and  we  denote  the  result  as  Fshft-input.  In  MATLAB  source  code,  this  may  be  written  as 
a  single  line  code:  fftshift(fft2(fftshift(5^-/npwt))).  Figures  4a  and  b  show  respectively 
the  input  excitation  function  and  the  result  after  applying  a  fftshift.  Figure  5b  shows  the 
result  after  applying  a  fft2  and  Figure  6b  shows  the  final  result  of  the  fast  Fourier 
transform  operation  after  applying  another  fftshift.  (Figure  6  is  shown  in  two- 
dimensional  perspective  so  that  negative  values  may  be  seen.)  The  absolute  value  is 
shown  in  Figure  7b.  (Because  these  graphs  are  cylindrically  symmetric  in  shape,  the  two- 
dimensional  perspective  here  will  illustrate  better  that  both  Figure  7a  and  b  are  equal.) 

Figures  5  a  and  6a  show  the  wrong  way  of  executing  a  fast  Fourier  transform  on 
an  input  excitation  function.  If  a  fft2  is  applied  directly  onto  (see  Figure  5a), 
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(a)  SHFT-INPUT  (b)  INPUT 


Y-axis  °  °  X-axis  Y-axis  ®  °  X-axis 


Figiire  4.  (a)  SHFT-INPUT  is  a  center  geometry  circular  excitation  and  (b)  INPUT  is  a 
comer  geometry  circular  excitation  obtained  by  applying  fftshift  to  the  center  geometry 

circular  excitation. 


(a)  FFT2(SHFT-INPUT)  (b)  FFT2(INPUT) 


fy-axis  ^  ®  fic-axis  fy-axis  ®  fx-axis 


Figure  5.  (a)  After  applying  fft2  on  Figure  4a  and  (b)  after  applying  fft2  on  Figure  4b. 
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(a)  FFTSHIFT(FFT2(SHFT-INPUT))  (b)  FFTSHIFT(FFT2(INPUT)) 


Figure  6.  (a)  After  applying  fTtshift  to  Figure  5a  and  (b)  after  applying  fTtshift  to  Figure 
5b.  These  graphs  are  shown  in  two-dimensional  perspective  so  that  negative  values  may 

be  seen. 


(a)  ABS(FFTSHIFT(FFT2(SHFT-1NPUT))) 
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(a)  ABS(FFTSHIFT(FFT2(INPUT))) 


Figure  7.  (a)  Absolute  value  of  Figure  6a  and  (b)  absolute  value  of  Figure  6b.  As  both 
graphs  are  cylindrically  symmetric  in  shape,  viewing  them  in  two-dimensional 
perspective  will  show  more  clearly  that  they  are  equal. 
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the  result  undergo  a  lateral  phase  shift  which  when  corrected  with  a  fftshift,  the  result  is 
shown  in  Figure  6a  and  its  absolute  value  in  Figure  7a.  Note  in  Figure  7  that  both 
methods  provide  a  similar  absolute  value  function  but,  as  shown  in  Figure  6,  the  Fourier 
transforms  are  different.  The  incorrect  method  gives  a  spiky  transform  and  has  wrong 
phase  information  required  for  our  computer  simulation. 

Beside  the  correct  order  of  application  of  these  two  fft  functions,  we  would  also 
like  to  highlight  another  very  important  fact  pertaining  to  their  speed  of  computation.  The 
MATLAB  User  Guide  [Ref  7]  points  out  that  when  the  row  and  column  dimensions  of 
the  matrix  are  power  of  two,  a  high-speed  radix-two  fft  algorithm  is  used.  When  the 
dimensions  are  not  other  than  a  power  of  two,  a  non-power-of-two  algorithm  finds  the 
prime  factors  of  the  dimensions  and  computes  the  mixed-radix  discrete  Fourier  transform. 
This  latter  process  can  be  quite  time  consuming,  particularly  as  the  size  of  the  matrices 
becomes  larger.  For  this  reason,  a  decision  was  made  to  work  with  NxN  matrices,  where 
Nisa  power  of  two. 

Now  that  we  have  discussed  the  theories  involved  in  this  thesis,  next  we  shall 
show  how  we  simulate  our  propagation  model  in  MATLAB. 
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III.  MATLAB  SIMULATION 


This  chapter  describes  the  simulation  programs  written  in  MATLAB.  Simulation 
is  used  here  to  refer  to  the  modeling  of  Equations  22  and  23  of  Chapter  II  in  MATLAB 
source  codes  and  the  animation  of  the  behavior  of  the  propagation  spatial  filter  and  the 
output  field.  Section  A  discusses  the  program  structure  adopted  for  our  simulations 
programs  and  section  B  explains  critical  algorithms  in  each  program  module.  No  in-depth 
knowledge  of  MATLAB  is  assumed  and  the  discussion  of  the  program  will  be  as 
functional  as  possible.  All  MATLAB  source  codes  can  be  found  in  Appendixes  A  to  E. 

A.  PROGRAM  STRUCTURE 

In  an  effort  to  shorten  simulation  time,  a  modular  program  structure  has  been 
selected.  The  objective  is  to  separate  the  time-consmning  and  repetitive  calculation 
algorithms  into  separate  independent  modules  from  the  main  program.  Most  often,  these 
modules  will  only  be  executed  once  and  their  results  are  stored  into  data  files  to  be 
recalled  later  for  use  by  other  program  modules  during  the  simulation  process  or  for 
generating  three-dimensional  graphs  for  print  out. 

In  general,  we  may  characterize  our  programs  into  five  main  functional  types: 

( 1 )  To  create  input  excitation  field  distribution,  u( x,  y,0,  t)  =  s(x,y,0)5(t)  of 
Equation  1.  We  have  here  the  m-files:  CIRCLE.m,  SQUARE.m, 
CRCGAUS.m  and  CRCBESS.m.  These  generate  input  excitations  with 
circular,  square,  circularly  truncated  Gaussian  and  circularly  truncated 
Bessel  field  distributions,  respectively.  These  programs  are  written  as 
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function  m-files  with  two  or  three  required  input  arguments.  They  can  be 
executed  independently  by  just  calling  the  fimction  name  and  providing  it 
with  the  required  input  arguments  at  the  MATLAB  command  window. 
For  example,  CIRCLE(d,N)  will  create  a  circular  input  excitation  field 
distribution  with  a  diameter  of  d  units  based  on  a  square  base  of  N  units 
size.  In  addition,  these  function  m-files  may  also  be  executed  as  an 
embedded  function  in  a  script  m-file.  In  our  program  structure,  we  utilize 
these  function  m-files  in  both  ways,  which  we  will  elaborate  in  the  later 
sections. 

(2)  To  create  the  propagation  spatial  filter,  c(f^,fy,z,t)  of  Equation  23 .  This 

is  done  by  the  m-file,  lOPTFIL.m  (which  stand  for  Improved  OPTical 
FILter;  the  prefix,  "Improved"  is  added  to  differentiate  this  m-file  from  a 
previous  work  on  an  m-file  in  [Ref  3]).  lOPTFIL.m  is  written  as  a 
function  m-file,  whieh  generates  data  required  for  our  simulation  program. 
The  data  generated  by  lOPTFIL.m  is  stored  in  two  MATLAB  database 
files  named  as  OPTVAR.mat  and  PJlNxn.mat  (".mat"  is  a  file  extension 
used  by  MATLAB  database  files).  OPTVAR  stand  for  OPTical 
VARiables  and  it  stores  all  the  initialized  parameters  required  for 
subsequent  programs  computations.  PJlNxn  is  an  acronym  comprising  of 
P  which  stands  for  Propagation,  J1  for  the  Bessel  function  of  the  first  kind 
contained  in  the  filter  function,  N  for  the  square  matrix  size,  N,  used  to 
store  the  filter  function  data  and  n  for  the  integer  number  ranging  from  1 
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to  61  representing  the  time  slice  when  the  propagation  spatial  filter 
fimction  is  computed. 

(3)  To  compute  the  temporal  spatial  field  distribution  at  the  output  plane, 
p(x,y,z,t).  This  is  done  by  lOPTPROP.m  (which  stand  for  Improved 
OPTical  PROPagation  for  the  same  reason  stated  in  the  above  paragraph). 
lOPTPROP  is  also  written  as  a  function  m-file  and  it  generates  data 
required  for  our  animation  programs.  The  main  function  of  lOPTPROP.m 
is  to  compute  p(x,y,z,t)  by  taking  the  inverse  Fourier  transform  of  the 
product  of  's(f^,fy)  and  c(f^,fy,z,t).  The  results  are  stored  in  two 

forms  in  two  separate  MATLAB  database  files,  OPTABS.mat  and 
OPTOUT.mat.  OPTABS.mat  contains  the  output  field  intensity,  which  we 
will  use  later  in  our  animation  programs  to  simulate  the  image  on  the 
output  plane.  OPTOUT.mat  contains  the  data  required  to  plot  a  three- 
dimensional  graphical  representation  of  the  total  output  field. 

(4)  To  plot  two-  and  three-dimensional  graphical  representations  of  the  input 
excitation  distribution,  the  propagation  spatial  filter  function  and  the 
output  field  distribution  in  both  temporal  and  spatial  frequency  domains. 
This  is  done  by  PLOTFILTER.m  and  PLOTFIELD.m  fi'om  data  generated 
by  lOPTFIL.m  and  IOPTPROP.m.  Both  of  these  programs  are  written  as 
script  m-files.  The  graphs  generated  by  these  two  programs  allow  us  to 
view  the  input  and  output  field  distributions  as  well  as  the  filter  function 
behavior  at  different  time  slices,  in  different  viewing  perspective  and  with 
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different  magnification  factors.  Hence,  they  provide  us  a  very  useful 
means  to  analyze  our  optical  propagation  model  at  different  stages  of  time 
and  space. 

(5)  To  animate  our  optical  propagation  model.  This  is  done  by  ANIMATE  1  .m 
ANIMATE2.m  and  ANIMATES  .m.  These  programs  are  written  as  script 
m-files.  Their  purposes  are  to  animate  the  behavior  of  our  propagation 
spatial  filter,  the  output  field  distribution,  the  total  output  field  distribution 
and  the  image  (or  field  intensity)  over  the  entire  simulation  time.  The  three 
animation  m-files  provide  similar  types  of  information  but  in  different 
formats.  This  is  to  cater  to  different  purposes  which  we  will  elaborate 
further  in  the  subsequent  sections. 

Figure  8  is  a  flow  chart  that  shows  our  program  structure  as  described  above  and 
from  this  figure,  we  may  see  the  inter-links  between  the  various  modules. 

B.  PROGRAM  DESCRIPTION 

The  following  sub-sections  explain  in  detail  the  critical  algorithms  in  each  of  the 
five  functional  file  types  discussed  above. 

1.  Input  Excitation  Field  Distribution  Program  Module 

As  mentioned  before,  the  m-files  required  to  generate  the  input  excitation  field 
distribution  comprise  of  CRCLE.m,  SQUARE.m,  CRCGAUS.m  and  CRCBESS.m. 

In  this  section,  we  would  like  to  highlight  the  reason  for  selecting  these  four 
specific  shapes  for  our  input  excitation.  They  were  selected  because:  (1)  they  represent 
real  input  excitation  sources  that  can  be  easily  generated  in  the  optical  laboratory  and  (2) 
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►I  3D  graph  filter  (PLOTFILTER.m)|^ 


Figure  8.  Programs  structure  and  program  flow. 


their  geometrically  symmetric  properties  help  to  simplify  our  MATLAB  algorithm  used 
to  generate  them.  The  latter  means  that,  since  the  shape  of  the  input  excitation  is 
symmetrical,  we  are  able  to  reproduce  the  whole  excitation  field  by  simply  generating  a 
quarter-shape  of  the  field  and  duplicating  it  three  times  to  create  the  whole  field. 
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However,  as  we  shall  see  later,  the  implementation  in  MATLAB  is  not  as  straight 
forward  as  discussed  here. 

For  reasons  already  explained  in  Chapter  I,  we  would  like  to  fix  the  size  of  the 
input  and  output  planes.  Therefore  in  our  program,  we  represent  these  planes  with  a 
square  base  matrix  of  size  N,  where  N  must  be  power  of  two  as  explained  at  the  end  of 
Chapter  II. 

In  MATLAB,  matrix  indices  begin  with  1  rather  than  0  (i.e.,  the  upper  left  entry 
being  row  1,  column  1  and  not  row  0,  column  0  as  an  origin  would  require).  This  means  a 
matrix  of  dimension  NxN  will  have  N  points  and  N-\  segments  in  each  row  and  column. 
This  also  means  that  the  center  of  symmetry  of  the  array  which  we  denote  as  NO,  would 
be  at  the  number  (N+l)/2  row  and  the  {N+\)/2  column.  However,  in  our  case,  we  require 
that  N  must  be  power  of  two  and  we  have  chosen  the  number  64  for  preliminary 
simulation  programs  (128  later  was  used  to  achieve  better  resolution).  This  means  we  will 
not  be  able  to  find  an  associate  center  of  symmetry  since  for  N=64,  (N+\)/2  will  not  be  an 
integer  number.  Therefore,  in  our  program,  we  had  to  arbitrary  choose  NO  to  be  as  near 
to  the  actual  center  as  possible  and  position  (33,33)  was  the  best  choice. 

Figure  9  depicts  a  NxN  equal  a  64x64-array  base  situated  on  the  x,  y  plane  divided 
into  four  quadrants.  To  generate  the  whole  input  excitation,  first  we  generate  the  field  to 
fill  the  smallest  quadrant,  which  is  quadrant  IV.  Then  by  flipping  up,  we  create  the  field 
in  quadrant  II  and  then,  by  flipping  both  quadrants  II  and  IV  to  the  left,  we  create  the 
field  in  quadrants  I  and  III,  thus  completing  the  field  on  the  entire  input  plane.  In 
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MATLAB,  the  flipping  of  the  field  can  be  achieved  by  using  the  flipud  and  fliplr 
commands  as  illustrated  in  the  source  code  of  Appendix  A. 
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Figure  9.  Base  array  configuration.  The  small  arrows  show  the  direction  of  flipping. 
Figure  10  shows  the  graphical  plot  of  the  four  input  excitation  field  distributions 
generated  by  the  above  program  module. 
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(a)  Circular  excitation 


Y-axis  ^  X-axis 


(b)  Square  excitation 


Y-axis  ®  ®  X-axis 


Figure  10.  Input  excitation  field  distribution  with  A^=  64:  (a)  circular  field  distribution 
with  if  =  25,  (b)  square  field  distribution  with  w  =  25,  (c)  circularly  truncated  Gaussian 
field  distribution  with  d  =25  and  a  =  1,  and  (d)  circularly  truncated  Bessel  field 
distribution  with  if  =  25  and  sigma  =12. 


28 


2.  Propagation  Spatial  Filter  Program  Module 

The  propagation  spatial  filter,  c(f^,fy,z,t)  is  computed  by  the  m-file 

lOPTFIL.m.  In  addition,  lOPTFIL.m  is  also  responsible  for  initializing  and  storing  all  the 
defining  parameters  into  a  data  file  named  OPTVAR.mat  that  is  required  for  other 
program  modules.  First,  we  will  go  through  all  defining  parameters  store  in 
OPTVAR.mat  and  then  we  will  explain  how  c(f^,fy,z,t)  is  computed  in  lOPTFIL.m. 

A  defining  parameter  is  a  parameter  that  delineates  an  aspect  of  the  basic  setup 
which  all  the  remaining  parameters  or  variables  depend.  Table  1  shows  all  the  defining 
parameters  used  in  this  thesis  and  their  assigned  values  used  for  our  simulation  model. 
The  first  parameter,  N,  sets  the  dimension  of  the  square  base  array  giving  the  number  of 
spatial  sample  points.  The  next  parameter,  NO,  defines  the  center  of  this  square  base 
array.  M  is  the  number  of  time  samples  or  time  slices,  which  we  use  to  observe  the  filter 
behavior  over  time.  Step  represents  the  number  of  leading  zeros  in  the  NxM  ontpvX  array. 
Step  is  required  to  simulate  the  Heaviside  step  function  that  we  see  in  Equation  21. 
Timejnax  is  the  maximum  propagation  time  that  we  have  fixed  for  our  simulation  model, 
z  is  the  distance  between  the  input  and  output  planes  and  c  is  the  speed  of  light.  We  have 
also  parameter,  rho,  which  represents  the  maximum  spatial  radius  of  the  filter  function. 

At  the  beginning  of  lOPTFIL.m,  we  initialize  all  the  above  defining  parameters  to 
the  assigned  values  for  our  simulation  model.  Then  we  used  some  of  these  parameters  to 
generate  two  important  matrices,  time  and  row,  which  are  required  for  the  program 
module,  IOPTPROP.m  for  the  computation  of  p(x,y,z,t)anA.  lOPTFIL.m  itself  for  the 
computation  of  c(f,^,fy,z,t)  respectively. 
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PARAMETER 

VALUE 

DEFINITION 

N 

64 

size  of  square  matrix 

NO 

33 

assigned  center  of  square  matrix 

M 

64 

total  number  of  time  slices 

Step 

3 

time  increment  prior  to  z/c 

timejnax 

0.95e-9  ns 

maximum  observation  time 

z 

100  mm 

distance  between  input  and  output  planes 

c 

3e8  m/s 

velocity  of  propagation 

rho 

200  mm 

spatial  radius 

Table  1.  Defining  parameters  and  their  assigned  values  for  our  propagation  model. 


The  matrix,  time,  represents  the  time  base  that  we  used  to  observe  the  filter 
function  and  the  output  field,  time  is  generated  by  the  MATLAB  built-in  function, 
linspace  with  z/c,  timejnax  and  M-Step  as  input  arguments.  Basically,  linspace  divides 
the  time  period  from  z/c  to  timejnax  into  M-Step  points  which  represent  61  time  slices, 
each  of  10  picoseconds  when  M  =  64  (if  M  =  128,  125  time  slices  each  of  5  picoseconds 
were  used).  These  ultrashort  time  slices  are  required  in  order  to  capture  the  fast  rate  of 
change  of  the  field  distribution  fi’om  the  pulsed  input  excitation.  With  these  ultrashort 
time  slices,  we  are  also  able  to  observe  our  propagation  model  in  slow  motion  in  our 
animation  program  modules,  which  no  existing  optical  measuring  equipment  is  capable 
of  doing. 

Next,  lOPTFIL.m  computes  the  filter  spatial  radius  matrix,  row  (which  is  actually 
the  parameter,  p,  of  Equation  23)  with  three  separate  steps.  First,  we  divide  the  value  of 
rho  into  NO-\  linear  spaces  with  the  linspace  command.  These  linear  spaces  are  then 
stored  in  the  matrix,  rhojn,  and  represent  the  radial  discrete  points  between  the  center  to 
the  side  of  the  NxN  matrix  space.  Then  the  cartesian  equivalent  of  rho  (with  coordinates 
label  as  rhox  and  rhoy)  is  generated  by  using  the  meshgrid  command  on  rhojn.  Next, 
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using  the  coordinates  rhox  and  rhoy  and  applying  the  Pythagorean  equation,  we  compute 
the  radial  distances  from  the  center  (i.e.,  NO)  to  any  points  on  one  of  the  four  quadrants 
of  the  NxN  matrix  space  (depicted  in  Figure  9).  These  radial  distances  are  then  stored  into 
a  NOxNO  matrix,  call  row.  Note  again  that  row  here  contains  only  the  radial  distances  for 
just  one  quadrant  of  the  NxN  matrix  space.  To  compute  the  radial  distances  for  the  whole 
NxN  matrix  space,  we  use  a  similar  algorithm  as  illustrated  in  Figure  9.  However,  as  we 
shall  see  later  that  we  do  not  apply  this  algorithm  straight  away  onto  row,  but  as  part  of 
the  computation  of  c( f^,fy,z,t) .  Once  we  have  computed  time  and  row,  we  store  the 
parameters,  N,  NO,  M,  Step,  c,  z  and  time  into  the  data  file,  OPTVAR.mat  with  the  save 
command  and  then  proceed  to  compute  c(f^,fy,z,t). 


As  given  by  Equation  23,  c(f,.,fy,z,t)  are  defined  over  three  different  time 


regions: 

(1)  t<  z/c.  This  represents  the  time  when  the  laser  pulse  has  not  yet  reached 
the  output  plane  and  hence  the  field  at  the  output  plane  is  zero.  Therefore,  we  will 
not  even  consider  this  time  region  in  our  simulation  model. 

(2)  t  =  z/c.  This  represents  the  time  when  the  laser  pulse  has  first  reached  the 
output  plane  and  hence  we  expect  a  sudden  jump  in  the  field  amplitude  given  by 
(as  in  Equation  23) 


c(f,.fy,z,t  =  zlc)  =  -zp  + 


2  _  2z/o(Q)  _  _  2  --  2 


2z 


=  -zp  = 


(24) 


In  our  simulation  program,  this  represents  time  slice  1. 
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use  Equation  25.  Note  that  when  we  compute  c(f^,fy,z,t)  by  using  Equation  24  or  25, 

we  are  only  computing  the  field  for  one  of  the  four  quadrants  of  the  l^xN  matrix  space.  To 
compute  the  whole  field  of  the  NxN  matrix  space,  we  adopt  the  same  algorithm  as 
depicted  in  Figure  9.  For  each  time  slice  when  c( f^,fy,z,t)  is  computed,  we  store  its 

value  into  a  matrix  named  as  PROP(m),  where  m  is  the  m*  time  slice  and  then  store  this 
matrix  into  a  data  files  PJlNxn.mat.  We  have  also  incorporated  a  "movie  play"  feature 
into  this  program  loop  to  allow  us  to  observe  the  changing  behavior  of  the  filter  function 
at  different  time  slices  with  the  command,  moviein,  getframe  and  movie.  The  source 
code  of  lOPTFIL.m  can  be  found  in  Appendix  B. 
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3.  Temporal  Spatial  Field  Distribution  Program  Module 

The  temporal  spatial  field  distribution,  p(x,y,z,t),  is  computed  by 
lOPTPROP.m.  IOPTPROP.m  also  uses  p(x,y,z,t )to  compute  the  field  intensity  on  the 
output  plane  by  taking  the  square  of  its  absolute  value  and  we  have  used  this  to  simulate 
the  image  that  would  be  seen  from  behind  the  output  plane. 

We  begin  the  computation  of  p(x,y,z,t)  by  first  loading  all  the  defining 
parameters  from  the  file  OPTVAR.mat  with  the  load  command.  Then  we  allow  the 
program  user  to  select  any  one  of  the  four  input  excitations:  Circle,  Square,  Gaussian  and 
Bessel.  The  Circle  and  Square  are  equal  amplitude  sources  having  the  shape  of  a  circle 
and  square,  respectively.  The  Gaussian  and  Bessel  inputs  are  circularly  truncated 
functions  that  have  spatially  varying  amplitudes  across  the  circular  face  of  the  source. 
After  the  program  user  has  selected  the  input  excitation,  the  user  is  asked  to  input  the 
diameter,  d,  of  the  truncating  circle  or  the  width,  w,  of  the  square  in  the  case  of  the 
Square  function.  For  the  case  of  the  Gaussian  and  Bessel  inputs,  the  user  is  further 
requested  to  input  the  standard  deviation,  sigma,  or  a  scaling  factor,  a,  respectively.  The 
purpose  of  this  part  of  the  program  is  to  not  only  allow  the  user  to  select  one  of  four 
choices  of  input  excitations  but  also  to  give  the  user  the  flexibility  to  vary  the  cross- 
sectional  size  of  the  excitation.  This  feature  of  the  program  has  expanded  the  user  choice 
to  analyze  and  study  varieties  of  input  excitations. 

We  shall  denote  the  selected  input  excitation  as  shft-input.  From  shft-input,  input 
is  created  by  shifting  shft  input  with  the  fftshift  command  to  a  comer  geometry.  Then 
with  the  fft2  command,  we  create  its  two-dimensional  spatial  Fourier  transform,  F-input. 
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As  explained  in  Chapter  II,  the  fftshift  operation  is  necessary  before  the  spatial  Fourier 
transform  operation  to  obtain  the  correct  phase  relationship  in  the  transform  operation. 
With  a  shift  back  to  the  center  geometry  with  another  fftshift  command,  the  angular 
spectrum  of  the  source  s(f^,fy),  called  Fshftjnput  in  the  program,  is  created.  The 
Bessel  ftmction  propagation  transfer  function,  c( Equation  23  must  now 
be  loaded  from  the  data  file,  P J 1  Nxn.mat.  The  product  of 's(f^,fy)  and  c(f^,f^,z,t)  is 
then  taken  to  find  p(f^,fy,z,t)  which  is  called  Fshft  output.  The  loading  and 

multiplication  process  is  repetitive  since  Fshftjnput  must  form  a  product  with  the  filter 
function  for  each  time  slice.  This  repetitive  multiplication  is  accomplished  with  a 
program  loop. 

To  find  the  desired  result  (i.e.,  p(x,y,z,t)),  the  two-dimensional  inverse  spatial 
Fourier  transform  (ifft2)  must  be  taken  for  the  product.  Before  this  can  be  done,  Fjoutput 
is  formed  by  shifting  Fshfi_output  from  the  centered  geometry  to  the  comer  geometry. 
Executing  the  inverse  transform  of  the  product  yields  output,  which  is  then  shifted  to  give 
shft_output.  The  array  shft_output  represents  the  output  at  the  time  slice  that  the  loop  is 
currently  computing.  (Note  that  shft_output  does  not  depict  the  optical  wave  or  the 
propagation  pattern  through  time;  it  only  depicts  the  optical  wave  at  a  specific  time). 

Because  of  the  cylindrical  symmetry  of  the  output  field,  to  produce  a  time  history 
of  the  desired  output  (which  we  refer  subsequently  as  the  "total  output"),  first  we  take  its 
absolute  value,  call  it  shft_outabs,  and  then  store  the  center  row  (row  NO)  of  shft_outabs 
into  another  array,  output _plot,  as  its  m  column  (where  m  is  the  loop  counter  which 
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relates  directly  to  the  time  slice  number).  For  example  ii  m  =  A,  then  column  4  of 
output  j)lot  represents  the  center  row  of  the  absolute  value  of  the  output  field  computed 
at  the  fourth  time  slice.  The  array  output _plot  is  therefore  of  size  NxM  and  we  store  this 
array  into  a  data  file  named  as  opdxM.mat  where  d  is  the  diameter  of  the  input  excitation 
and  Mis  the  number  of  time  slices. 

To  obtain  the  field  intensity,  the  square  of  shft_output  is  taken  and  this  is  called 
shftjntensity.  For  each  time  slice,  shftjntensity  is  stored  in  a  file  named  optabm.mat 
where  m  corresponds  to  the  /w*  time  slice  when  the  field  intensity  was  computed.  We  use 
optabm.mat  in  the  animation  program  modules  to  simulate  the  image  at  the  output  plane. 
The  source  code  for  lOPTPROP.m  can  be  found  in  Appendix  C. 

4.  Two-  and  Three-Dimensional  Graphical  Program  Modules 

There  are  two  program  modules  that  do  two-  and  three-dimensional  graphical 
plotting  of  the  filter  function  and  the  output  field  at  different  time  slices.  These  are  called 
PLOTFILTER.m  and  PLOTFIELD.m.  While  PLOTFILTER.m  makes  use  of  the  data 
from  data  files,  PJlNxn.mat,  to  plot  the  filter  function  at  different  time  slices  in  different 
perspectives,  PLOTFIELD.m  makes  use  of  the  data  from  data  files,  opdxM.mat  and 
optabm.mat,  to  plot  the  output  field  and  total  output  at  different  time  slices  in  different 
perspectives.  The  basic  commands  used  in  these  two  program  modules  are  load,  figure, 
mesh,  subplot,  axis,  grid,  set  and  view.  Of  all  these  commands,  view  is  one  of  the  most 
useful  ones  as  we  use  view  in  our  programs  very  frequently  to  plot  graphs  in  different 
perspectives.  Figures  11  and  12  show  examples  of  some  three-dimensional  graphical 
plots  generated  by  these  two  program  modules. 
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Time  slice  2:  Filter  response 
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Time  slice  20:  Filter  response 
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Figure  1 1 .  Three-dimensional  graphs  of  the  filter  fimction  (left)  and  output  field  (right)  at 
time  slices  1,2, 10  and  20.  Notice  that  as  the  time  slice  number  increases  the  amplitude 
of  the  field  decreases  but  the  field  radial  spreading  increases. 
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Time  slice  30:  Filter  response 
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Figure  12.  Three-dimensional  graphs  of  the  filter  fimction  (left)  and  output  field  (right)  at 
time  slices  30, 40,  50  and  61 .  Notice  that  as  the  time  slice  number  increases  the  amplitude 
of  the  field  decreases  but  the  field  radial  spreading  increases. 


The  graphs  on  the  left  represent  the  filter  function  and  those  on  the  right  represent  the 
output  field  in  different  time  slices.  The  source  code  of  PLOTFILTER.m  and 
PLOTFIELD.m  can  be  found  in  Appendix  D. 

5.  Animation  Program  Modules 

In  addition  to  the  usual  static  two-  and  three-dimensional  plots  generated  by  the 
graphical  program  modules  described  above,  the  animation  program  modules  go  one  step 
further  to  animate  the  changing  behavior  of  the  filter  as  well  as  the  output  fields  at  the 
output  planes  over  the  entire  simulation  time.  We  have  adopted  the  "ffame-by-ffame 
capture  and  playback"  technique  to  create  our  animation.  For  example,  to  animate  the 
filter  function,  we  plot  the  filter  function  at  each  time  slice,  starting  from  time  slice  1  to 
61 .  For  each  time  slice  that  we  plot,  this  represents  a  frame  of  our  animated  movie  and,  if 
we  were  to  capture  and  playback  all  the  frames  in  sequence,  it  gives  the  effect  of  a 
moving  picture. 

Three  main  commands  are  used:  moviein,  getframe  and  movie.  To  create  an 
animation,  the  command  moviein  is  first  used  to  pre-allocate  enough  memory  space  to 
store  all  the  graphical  frames  which  comprises  of  the  61  PROP(m)  filter  function,  61 
shft_output  and  61  shft  intensity  matrices  computed  over  the  61  (M-Step)  time  slices. 
This  is  a  lot  of  memory  especially  when  we  are  plotting  three-dimensional  graphs. 
Therefore,  to  save  on  computer  memory  and  thus  to  increase  the  speed  of  animation,  we 
combine  all  three  matrices  into  a  single  graphical  plot  using  the  subplot  command  and 
capture  all  three  plots  into  a  single  frame.  In  this  way,  we  need  only  to  pre-allocate 
enough  memory  for  61  frames.  To  increase  the  speed  of  animation  further,  we  have  also 
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reduced  the  size  of  the  viewing  window,  which  further  reduces  the  computer  memory 
required  to  do  the  animation. 

We  use  the  getframe  command  to  take  a  snapshot  of  the  current  plot  for  use  in 
the  movie  playback  and  to  start  playback,  we  use  the  movie  command.  Three  animation 
program  modules  were  written  and  they  all  animate  the  filter  function,  the  output  field 
and  the  image  at  the  output  plane  but  in  different  format.  These  formats  are  shown  in 
Figures  13, 14  and  15. 

Now  that  we  have  explained  in  detail  all  the  program  modules,  we  shall  proceed 
to  the  next  chapter  on  numerical  simulation. 


Figure  1,3.  Animation  format  1.  The  graph  (a)  shows  the  filter  spatial  frequent  response, 
graph  (b)  shows  the  output  field  and  graph  (c)  shows  the  image  at  the  output  plane.  There 
is  a  small  window  just  below  graph  (c)  that  shows  the  time  slice. 
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(a)  Filter  spatial  frequency  response 


(b)  Image  field  intensity 
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Figure  14.  Animation  format  2.  Graph  (a)  shows  the  filter  spatial  frequent  response, 
graph  (b)  shows  the  output  field,  graph  (c)  shows  the  close-up  side  view  of  the  total 
output  and  graph  (d)  shows  the  image  at  the  output  plane.  There  is  a  small  window  just 

below  graph  (d)  that  shows  the  time  slice. 
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(a)  Filter  spatial  frequency  response 
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axis 


(b)  Image  field  intensity 


Y-axis  0  0  X-axis 


(c)  Field  distribution 


Time-slice 


8 

6 

4 

2 


X-axis 


time-slice  61 


Figure  15.  Animation  format  3.  Graph  (a)  shows  the  filter  spatial  frequent  response, 
graph  (b)  shows  the  output  field,  graph  (c)  shows  a  ten  times  magnified  view  of  the  total 
output  and  graph  (d)  shows  the  image  at  the  output  plane.  There  is  a  small  window  just 

below  graph  (d)  that  shows  the  time  slice. 
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IV  NUMERICAL  SIMULATION 


This  chapter  presents  the  numerical  simulation  results  for  the  MATLAB  program 
modules  described  in  Chapter  III.  Section  A  shows  and  analyzes  the  propagation  spatial 
filter  fimction  over  different  time  slices.  Section  B  steps  through  the  process  for 
generating  the  results  for  a  circular  field  input  excitation.  An  analysis  is  also  done  here  on 
the  effect  of  varying  the  diameter  of  the  circular  field.  Finally,  the  simulation  results  for 
the  square,  circularly  truncated  Gaussian  and  circularly  truncated  Bessel  input  excitations 
are  presented  and  compared  with  the  circular  input  field  in  term  of  the  field  amplitude 
and  rate  of  field  radial  spreading.  We  also  discuss  the  formation  of  constructive 
interference. 

A.  PROPAGATION  SPATIAL  FILTER  FUNCTION 

IOPTFIL.m  was  first  executed  with  W=  64  and  M=  64  and  Figure  16  provides 
the  results  of  the  spatial  filter  function  for  time  slices  1,  2,  30  and  61.  (These  time  slices 
were  arbitrary  chosen  to  illustrate  the  dynamic  changes  of  the  filter  over  the  entire 
simulation  time.)  The  waveform  at  time  slice  1  represent  the  filter  function  computed  at  t 
=  z/c  and  the  waveforms  at  time  slices  2,  30  and  61  represent  the  filter  function  computed 
when  t  >  z/c. 

Notice  that  as  the  time  slice  increases,  the  filter  fimction  becomes  very  coarse  and 
spiky  (notably  at  the  higher  time  slices  such  as  30  and  61  as  shown).  In  order  to  smooth 
these  waveforms,  we  need  a  higher  number  of  spatial  sampling  points  and  hence  we  set  N 
=  128.  With  this  new  value  of  N,  we  are  able  to  observe  finer  details  of  the  changes 
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occurring  in  the  filter  fimction  as  time  progresses.  Furthermore,  we  have  also  set  M=  128 
so  that  we  may  observe  a  more  progressive  change  in  the  filter  function. 


Figure  16.  Propagation  spatial  filter  function  with  64  and  M=  64:  (a)  time  slice  1,  (b) 
time  slice  2,  (c)  time  slice  30  and  (d)  time  slice  61.  Note  that  time  slice  1  occurs  at  t  =  z/c 
and  time  slices  2,  30  &  61  occur  at  t  >  z/c.  Notice  that  these  waveforms  look  very  coarse 

and  spiky. 
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The  resultant  filter  fiinction  obtained  with  the  new  set  of  values  for  N  and  M  are 


shown  in  Figures  17  and  18.  Figure  17  shows  the  filter  fiinction  at  time  slices  1,  2,  4  and 
8.  Note  that  both  the  waveforms  at  time  slice  1  of  Figure  16  and  17  are  similar  as  they 
both  represent  the  filter  function  computed  at  t  =  z/c. 


However,  the  waveforms  at  time  slice  2  of  Figure  16  and  17  are  not  similar  because  the 
time  slice  increment  for  both  waveforms  are  not  the  same;  recall  from  Chapter  III  that 
when  M=  64,  each  time  slice  increment  is  10  picoseconds  whereas  for  M=  128,  each 
time  slice  increment  is  only  5  picoseconds.  Figure  1 8a  and  b  show  the  propagation  spatial 
filter  function  at  time  slices  60  and  125  respectively. 


(a)  Time  slice  60  (b)  Time  slice  125 
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(c)  Time  slice  60 
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(d)  Time  slice  125 


Figure  18.  Propagation  spatial  filter  function  with  7V^=  128  and  M=  128:  (a)  time  slice  60, 
(b)  time  slice  125,  (c)  filter  function  cross-section  view  at  time  slice  60  and  (d)  filter 
function  cross-section  view  at  time  slice  125.  Notice  that  more  peaks  are  formed  at  higher 

time  slice  number. 
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Notice  that  as  the  time  slice  number  increases,  the  filter  function  forms  more  peaks  and 
collapses  on  itself.  This  phenomenon  can  be  observed  by  looking  at  the  cross-sectional 
vie\v  of  the  filter  function  as  shown  in  Figure  18c  and  d  for  time  slices  60  and  125 
respectively. 

B.  OUTPUT  FIELD  DISTRIBUTION 

In  this  section,  we  shall  first  step  through  the  process  for  computing  the  output 
field  for  a  circular  field  input  excitation.  We  shall  present  results  obtained  for  a  circular 
source  with  diameter,  d=  25.  Another  set  of  results  was  also  generated  for  a  circular 
soxirce  vvith  d=A9to  compare  the  effect  of  an  increased  source  cross-sectional  area.  Then 
we  shall  present  the  simulated  results  for  the  square,  circularly  truncated  Gaussian  and 
circularly  truncated  Bessel  input  field  excitations. 

Before  we  proceed  further,  we  wish  to  highlight  that  the  diameter  for  all  the 
circular  soxirces,  d,  and  the  width  for  the  square  source,  w,  are  given  in  term  of  the 
number  of  spatial  points  on  the  input  planes  which  is  made  up  of  a  NxN  array.  To  convert 
d  into  unit  of  centimeters  (given  by  D),  we  adopt  the  following  conversion  equations: 

D[cm]  =  100  *  d  *  Ak  ,  (26) 

and 


Ax  = 


1 

2  *  P„,a. 


1 

2  *  200 


^  0.0025 [m] 


(27) 


where  Ax  represents  the  lateral  displacement  between  each  discrete  points  on  the  NxN 
array  and  =  200  cycle/meter  is  the  spatial  radius  that  we  have  selected  for  our 
propagation  model.  For  example,  \fd=  25,  this  represents  a  physical  aperture  diameter  of 
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6.25  centimeters  on  the  input  plane.  In  what  follows,  we  shall  indicate  the  diameter,  D  in 
units  of  centimeter  in  parenthesis. 

1.  Circular  Field  Input  Excitation  with  Small  Aperture 

Figures  19  to  22  show  the  simulation  results  obtained  for  a  circular  field  input 
excitation  with  d- 25  (6.25  centimeter).  Figure  19a  shows  the  circular  input  excitation 
and  Figure  19b,  c  and  d  illustrate  the  three  steps  required  to  produce  its  two-dimensional 
spatial  Fourier  transform.  Note  that  Figure  19d  is  's(  f^,fy,0)  of  Equation  22  and  if  we 

multiply  this  with  c( f^,fy,z,t)  and  take  the  two-dimensional  inverse  spatial  Fourier 

transform  of  the  product,  we  produce  the  field  at  the  output  plane,  p(  x,y,z,t) . 

For  this  simulation  run,  we  use  M  =  128.  This  implies  that  we  have  125  time 
slices  and  the  simulation  run  have  computed  the  filter  function  125  times.  Unfortunately, 
we  are  unable  to  show  all  the  125  plots  of  the  filter  function  and,  therefore,  only  four  of 
these  plots  are  selected  to  illustrate  the  computation  of  the  output  field.  Figure  20  shows 
the  filter  function  computed  at  time  slices  1,  50,  100  and  125.  Notice  that  the  filter 
function  forms  more  peaks  and  collapses  on  itself  as  the  time  slice  number  increases.  If 
we  multiply  these  plots  with  that  of  Figure  19d  and  take  the  two-dimensional  spatial 
Fourier  transform  of  their  product  one  at  a  time,  we  produce  the  output  field  distribution 
as  illustrated  on  the  left-hand  side  plots  of  Figure  21.  Notice  that  the  amplitude  of  the 
output  field  decreases  as  expected  and  the  field  radial  size  increases  as  the  time  slice 
number  increases.  The  plots  on  the  right-hand  side  of  Figure  21  show  the  image  formed 
on  the  output  plane.  From  these  two-dimensional  plots,  the  field  radial  spread  with  time 
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becomes  more  obvious  and  this  has  clearly  demonstrated  the  effect  of  wave  diffraction 
over  time. 


(a)SHFT-INPUT  (b)  INPUT 


Y-axis  ^  ®  X-axis  Y-axis  ^  X-axis 


Figure  19.  Fourier  transform  of  an  impulse  plane  wave  illuminating  a  circular  aperture 
with  <:/=  25  (6.25  cm):  (a)  circular  input  excitation  field,  shft-input,  (b)  after  applying  a 
fftshift  on  shft-input  to  produce  input,  (c)  after  applying  fft2  on  input  to  produce  F-input 
and  (d)  after  applying  a  fftshift  on  F-input  to  produce  the  two-dimensional  spatial 
Fourier  transform,  Fshft-input.  Note  Fshft-input  represents  !( f,.,fy,0) . 
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To  present  the  output  field  distribution  on  the  output  plane  over  time  requires  a 
four-dimensional  plot  in  x-,  y-,  amplitude-  and  time-spaces.  However,  this  is  beyond 
MATLAB's  graphics  capability.  Nevertheless,  we  can  still  use  MATLAB's  three- 
dimensional  graphics  capability  to  show  this  four-dimensional  information  on  a  three- 
dimensional  plot.  To  do  this,  we  utilize  the  fact  that  the  output  field  distribution  is 
cylindrically  symmetric  and  therefore  a  center  cross-section  of  the  field  distribution  is 
sufficient  to  describe  the  entire  field  on  the  whole  output  plane.  If  we  line  up  all  the 
center  cross-section  of  the  output  field  for  all  the  125  time  slices,  we  are  able  to  produce  a 
time  history  of  the  output  field  and  we  denote  this  as  the  total  output. 

Figure  22a  shows  the  total  output.  As  expected,  we  observe  a  sudden  peak  at  time 
slice  1  where  z/c.  and  the  field  amplitude  is  given  by  Equation  24.  From  Figure  22b 
which  shows  a  ten  times  magnified  version  of  the  total  output,  we  observe  that  as  the 
time  slice  number  increases,  the  amplitude  decreases  and  eventually  goes  to  zero  in 
accordance  to  Equation  25.  Also,  we  can  witness  a  phenomenon  called  constructive 
interference  where  the  two  inboard  tails  meet  somewhere  at  time  slices  82.  This  manifest 
itself  as  an  unexpected  amplitude  peak  as  shown  in  the  close-up  cross-section  view  of  the 
total  output  in  Figure  22c.  From  Figure  22d,  which  shows  a  close-up  fi’ont  view  of  the 
total  output,  we  observe  the  increase  in  radial  spread  of  the  field  as  the  time  slice  number 
increases. 

2.  Circular  Field  Input  Excitation  with  Large  Aperture 

Figure  23  to  25  show  the  simulation  results  obtained  for  a  circular  field  input 
excitation  with  c?=  49  (12.25  centimeter).  With  a  d  value  of  almost  double  the  previous 
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(a)  SHFT-INPUT  (b)  FSHFT-INPUT 


Figure  23.  (a)  Circular  input  excitation  with  d  =  A9  (12.25  cm)  and  (b)  two-dimensional 

spatial  Fourier  transform. 

Figures  24a  and  b  show  the  two-dimensional  cross-section  view  of  the  two-dimensional 
spatial  Fourier  transform  for  d  =  25  and  d  =  49,  respectively.  From  this  figure,  we  may 
also  observe  that,  for  d  =  49,  in  addition  to  having  a  higher  peak  value,  the  input  also  has 
a  slimmer  spatial  Fourier  transform.  This  is  analogous  to  a  broad  time-base  signal  having 
a  narrower  spectral  response  than  a  narrow  time-base  signal. 


(a)  Small  circular  FSHFT-INPUT  (b)  Large  circular  FSHFT-INPUT 


Figure  24.  Two-dimensional  spatial  Fourier  transform  for  circular  input  excitation  for  (a) 

d=25  (6.25  cm)  and  (b)  =  49  (12.25  cm). 
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Figure  25  shows  the  total  output  for  d=  49  (12.25  centimeter)  in  different  perspective. 
From  these  plots,  we  can  observe  similar  phenomena  displayed  in  Figure  22  for  d  =  25. 


(a)  Total  output 


(b)  Magnified  view  of  total  output 


(c)  Close-up  cross-section  view  of  total  output 


(d)  Close-up  front  view  of  total  output 


space 


Figure  25.  Circular  field  input  excitation  with  49  (12.25  cm):  (a)  Total  output,  (b)  ten 
times  magnified  view  of  total  output,  (c)  close-up  cross-section  view  of  total  output  and 

(d)  close-up  front  view  of  total  output. 


However,  we  notice  that  in  this  case,  because  of  the  larger  cross-section  area  of  the 
wavefront,  the  radial  spreading  is  more  dispersed  and  that  the  strong  constructive 
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interference  that  we  observed  for  d=25  cannot  be  seen  here.  Nevertheless,  if  we  could 
extend  the  time  base  beyond  125  time  slices,  it  is  anticipated  that  the  two  inboard  tails 
would  still  meet  and  form  a  constructive  interference.  But  the  interference  amplitude 
would  probably  be  lower  because  the  expected  field  amplitude  beyond  time  slice  125  is 
also  lower. 

3.  Square  Field  Input  Excitation 

Figures  26  and  27  show  the  simulation  results  obtained  for  a  square  field  input 
excitation  with  w  =  25  (6.25  centimeter).  We  have  purposely  selected  w  =  25  so  that  we 
may  compare  this  set  of  results  with  that  obtained  for  the  circular  field  input  with  d  =  25. 
Notice  from  Figure  26  that  that,  because  a  square  input  of  width,  w  =  25,  has  a  larger 
cross-section  area  than  a  circle  input  with  d  =  25,  the  peak  amplitude  of  the  square's  two- 
dimensional  spatial  Fourier  transform  is  larger  (when  comparing  Figure  19d  with  Figure 
26b). 


(a)  SHFT-INPUT  (b)  FSHFT-INPUT 


Y-axis  ^  ®  X-axis  fy-axis  ®  b(-axis 


Figure  26.  (a)  Square  input  excitation  with  w  =  25  (6.25  cm)  and  (b)  two-dimensional 

spatial  Fourier  transform. 
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Figure  27  shows  the  total  output  for  the  square  input  excitation  with  w  =  25  (6.25 
centimeter)  in  different  perspective.  From  these  plots,  we  can  observe  similar  phenomena 
displayed  in  Figure  22  for  the  circular  input  excitation  with  d=  25. 


(a)  Total  output 


(b)  Magnified  view  of  total  output 


(c)  Close-up  cross-section  view  of  total  output 


(d)  Close-up  front  view  of  total  output 


space 


Figure  27.  Square  field  input  excitation  with  w  =  25  (6.25  cm):  (a)  Total  output,  (b)  ten 
times  magnified  view  of  total  output,  (c)  close-up  cross-section  view  of  total  output  and 

(d)  close-up  front  view  of  total  output. 
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time  slice 


For  this  case,  a  constructive  interference  has  occurred  at  time  slice  86  compared  to  82  for 
the  small  circular  input.  This  result  is  consistent  with  what  we  have  anticipated  for  the 
large  circular  input  that  constructive  interference  would  occur  at  further  time  slice  when 
the  wavefront  is  larger. 

4.  Circularly  Truncated  Gaussian  Field  Inpnt  Excitation 

Figures  28  and  29  show  the  simulation  results  obtained  for  a  circularly  truncated 
Gaussian  field  input  excitation  with  d=25  (6.25  centimeter)  and  sigma  =  12.  The  sigma 
factor  represents  the  standard  deviation  and  it  determines  the  width  of  the  full  Gaussian 
field.  Now,  because  a  truncated  Gaussian  field  has  a  wavefront  that  has  cross-section  area 
that  is  smaller  than  that  of  a  full  circular  input  field,  we  can  expect  to  see  a  lower  peak 
amplitude  for  its  two-dimensional  Fourier  transform  (when  comparing  Figure  19d  with 
Figure  28b). 


Figure  28.  (a)  Circularly  truncated  Gaussian  field  input  excitation  with  25  (6.25  cm) 
and  sigma  =12  and  (b)  Two-dimensional  spatial  Fourier  transform. 
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Figure  29  shows  the  total  output  for  the  circularly  truncated  Gaussian  field  input 
excitation  with  i/=  25  (6.25  centimeter)  and  sigma  =  12  in  different  perspective. 


(c)  Close-up  cross-section  view  of  total  output 


(d)  Close-up  front  view  of  total  output 


space 


Figure  29.  Circularly  truncated  Gaussian  field  input  excitation  with  d=25  (6.25  cm)  and 
sigma  =  12:  (a)  Total  output,  (b)  ten  times  magnified  view  of  total  output,  (c)  close-up 
cross-section  view  of  total  output  and  (d)  close-up  j&ont  view  of  total  output. 


From  these  plots,  again  we  can  observe  similar  phenomena  displayed  in  Figme  22  for  the 


circular  input  excitation  with  d  =  25.  Notice  that  in  this  case,  because  of  the  smaller 
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wavefront,  the  field  peak  amplitude  at  time  slice  1,  the  field  radial  spreading  and  the 
constructive  interference  amplitude  are  smaller  (when  comparing  Figure  22  with  Figure 
29). 

5.  Circularly  Truncated  Bessel  Field  Input  Excitation 

Figures  30  and  31  show  the  simulation  results  obtained  for  a  circularly  truncated 
Bessel  field  input  excitation  with  d=25  (6.25  centimeter)  and  o  =  1.  Factor  a  here  is  the 
width  scaling  factor  of  the  Bessel  fimction.  Figure  30a  and  b  show  the  circularly 
truncated  Bessel  field  input  excitation  and  its  two-dimensional  spatial  Fourier  transform 
respectively.  Notice  the  broad  spectral  response  in  Figure  30b  is  attributed  by  the  narrow 
waveform  of  the  input  field  (when  comparing  Figure  19  with  Figure  30). 

(a)  SHFT-INPUT  (b)  FSHFT-INPUT 


Y-axis  X-axis  fy-axis  °  °  fx-axis 

Figure  30.  (a)  Circularly  truncated  Bessel  field  input  excitation  with  25  (6.25  cm)  and 
a  =\2  and  (b)  two-dimensional  spatial  Fourier  transform. 

Figure  31  shows  the  total  output  for  the  circularly  truncated  Bessel  field  input  excitation 

with  d=2S  (6.25  centimeter)  and  a  =  1  in  different  perspective.  From  these  plots,  again 

we  can  observe  similar  phenomena  displayed  in  Figure  22  for  the  circular  input  excitation 


60 


with  d=  25.  Notice  also  that,  because  of  the  smaller  wavefront,  the  field  peak  amplitude 
at  time  slice  1  and  the  field  radial  spreading  are  smaller  (when  comparing  Figure  22  vdth 
Figure  31).  Moreover,  the  field  amplitude  is  decaying  so  fast  that  no  noticeable 
constructive  interference  is  created  at  the  point  where  the  two  inboard  tails  meet. 


(a)  Total  output 


(b)  Magnified  view  of  total  output 


(c)  Close-up  cross-section  view  of  total  output 


(d)  Close-up  front  view  of  total  output 


space 


Figure  31.  Circularly  truncated  Bessel  field  input  excitation  with  £?=  25  (6.25  cm)  and  a 
=  1 :  (a)  Total  output,  (b)  ten  times  magnified  view  of  total  output,  (c)  close-up  cross- 
section  view  of  total  output  and  (d)  close-up  front  view  of  total  output. 
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In  this  chapter,  we  have  done  a  detailed  analysis  of  the  characteristic  of  the 
propagation  spatial  filter  function.  We  have  also  presented  the  output  field  distribution 
for  the  four  different  shaped  input  fields  as  predicted  by  our  simulation  model.  With  this, 
we  have  come  to  the  concluding  chapter  of  our  thesis  research.  In  the  next  chapter,  we 
will  summarize  our  achievements  made  in  this  research  and  give  recommendations  for 
future  work. 
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V.  SUMMARY 


This  thesis  presented  a  MATLAB  implementation  of  a  Fourier  transform 
approach  to  model  and  predict  transient  optical  wave  propagation  through  tree-space. 
Linear  systems  theory  characterized  the  wave  propagation  model  in  terms  of  a  Green's 
function,  which  solves  the  wave  equation  and  satisfies  the  boundary  conditions  of  our 
propagation  model.  Fourier  transform  theory  simplified  the  mathematics  required  for  our 
computer  simulation. 

A  modular  programming  approach  was  adopted  for  our  MATLAB  program  to 
segregate  the  time-consuming  processes  from  the  less  time-consuming,  which  allow  the 
simulation  programs  to  run  more  efficiently  with  less  computer  memory.  User-interactive 
features  introduced  in  the  programs  allow  the  program  user  to  select  a  variety  of  input 
excitations  for  analysis.  Animation  programs  provided  visualization  of  the  changes  in  the 
filter  function  and  the  output  field  over  time.  Two-dimensional  plots  of  the  field  intensity 
were  presented  to  help  comprehend  the  image  formation  on  the  output  plane.  Detailed 
description  of  all  the  program  modules  were  given  and  Appendixes  A  to  E  contain  the 
source  codes.  Many  two-  and  three-dimensional  graphics  in  different  perspectives  were 
generated  to  demonstrate  the  program's  operation. 

We  computed  the  spatial  impulse  responses  for  a  circular,  square,  circularly 
tnmcated  Gaussian  and  circularly  truncated  Bessel  input  excitations.  Their  results  were 
compared  and  thoroughly  analyzed  to  identify  known  optical  phenomena  like  wave 
diffraction,  dispersion  and  constructive  interference. 
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Future  investigation  is  open  in  several  areas.  A  detailed  comparison  of  our 
simulation  model  should  be  made  with  existing  published  models.  The  physical 
interpretation  of  dit-zjc)  in  Equation  24  for  the  propagation  spatial  filter  should  be 
further  investigated,  especially  in  the  role  that  plays  our  time  plots  of  the  spatial  filter.  In 
this  thesis,  we  have  only  scratched  the  surface,  so  to  speak,  of  computer  simulation  of  the 
wide  ranging  techniques  of  optical  processing.  There  is  much  fascination  to  be  found  in 
aperture  design,  complex  filtering,  optical  computing,  etc.  Computer  graphics  allow  us  to 
visualize  more  than  the  eye  can  see. 
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APPENDIX  A.  SOURCE  CODE  FOR  INPUT  EXCITATIONS 
The  followings  are  the  MATLAB  source  code  for  the  four  input  excitations: 
circular,  square,  circularly  truncated  Gaussian  and  circularly  truncated  Bessel  field 
distributions. 

CIRCULAR  INPUT  EXCITATION 

function  Y  =crcle(d,N) 

%crcle.m:  Y=crcle(d,N) 

%Program  for  generating  uniform  circular  excitation  functions 
%d  is  the  diameter  of  the  circle.(ODD  integer) 

%N  is  the  width  of  the  square  base.  (EVEN  integer) 

%Example:  z  =  crcle(33,64); 

%  Adopted  from  [Ref  3] 

%Check  that  d  is  an  odd  integer 
if  rem(d,2)  <  0.1; 

error('The  diameter  of  the  circle  function  must  be  an  ODD  integer.'); 

else; 

end; 

%Check  that  N  is  an  even  integer 
if  rem(N,2)  ~=  0.0; 

error('The  width  of  the  square  base  must  be  an  EVEN  integer.'); 

else; 

end; 

NO  =  (N/2)+l ;  %NO  is  the  center  location 

r  =  d/2;  %r  is  the  radius 

Y  =  zeros(N); 
temp  =  zeros(NO-l); 
for  m=l:r+l; 

forn=l:r+l; 
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if  sqrt((m-l)^2  +  (n-l)^2)  <=  r; 
temp(m,n)=l ; 

end; 

end; 

end; 

%Generate  the  entire  N  x  N  input  function 
Y(NO:N,NO:N)  =  temp; 

Y(2;NO,NO:N)  =  flipud(temp); 
Y(2:NO,2:NO)  =  rot90(temp,2); 
Y(NO:N,2:NO)  =  fliplr(temp); 

%%%  End  of  program  %%% 


SQUARE  INPUT  EXCITATION 

function  Y  =  square(w,N) 

%square.m:  Y  =  square(w,N) 

%Program  for  generating  a  imiform  square  excitation  function. 

%w  is  the  width  of  the  table.  (ODD  inteqer) 

%N  is  the  width  of  the  square  base.  (EVEN  integer) 

%Example:  z  =  square(33,64); 

%  Adopted  from  [Ref.  3] 

%Check  that  w  is  an  odd  integer, 
if  rem(w,2)  <0.1; 

error('The  width  of  the  table  must  be  an  ODD  integer.'); 

else; 

end; 

%Check  that  N  is  an  even  integer 
if  rem(N,2)~=  0.0; 

error(’The  width  of  the  square  base  must  be  an  EVEN  integer.'); 
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else; 

end; 

NO  =  (N/2)+l ;  %NO  is  the  center  location 

wO  =  ceil(w/2);  %wO  is  the  mid-point  of  the  table 

Y  =  zeros(N); 
temp  =  zeros(NO-l); 
temp(l:wO,l:wO)=  ones(wO); 

Y(NO:N,NO:N)  =  temp; 

Y(2:NO,NO;N)  =  rot90(temp); 

Y(2:NO,2:NO)  =  rot90(temp,2); 

Y(NO:N,2:NO)  =  rot90(temp,3); 

%%%  End  of  program  %%% 


CIRCULARLY  TRUNCATED  GAUSSIAN  INPUT  EXCITATION 

function  Y  =  crcgaus(sigma,d,N) 

%crcgaus.m:  Y  =  crcgaus(sigma,d,N) 

%Program  for  generating  circular  Gaussian  functions. 

%sigma  is  the  standard  deviation  of  the  Gaussian  function. 

%d  is  the  diameter  of  circle.  (ODD  integer) 

%N  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 

%Example:  z  crcgaus(12,33,64); 

%  Adopted  from  [Ref.  3] 

mu=0;  %mu  is  the  mean  of  the  Gaussian  function 

%Check  that  d  is  an  odd  integer. 
ifrem(d,2)  <0.1; 

error('The  diameter  of  the  circle  function  must  be  an  ODD  integer.'); 

else; 

end; 
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%Check  that  N  is  an  even  integer, 
if  rem(N,2)  ~=  0.0; 

error(The  width  of  the  square  base  must  be  an  EVEN  integer.'); 

else; 

end; 

NO  =  (N/2)+l ;  %NO  is  the  center  location 
r  =  d/2;  %r  is  the  radius 

Y  =  zeros(N); 
temp  =  zeros(NO-l); 
form=  l:(d+l)/2; 

forn=  l:(d+l)/2; 

X  =  sqrt((m- 1  )^2  +  (n- 1  )''^2); 
ifx<=r; 

temp(m,n)  =  (l/(sqrt(2*pi)*sigma))*exp(-((x-mu)^2)/(2*(sigma^2))); 
end; 

end; 

end; 

Y(NO:N,NO:N)  =  temp; 

Y(2:NO,NO:N)  =  flipud(temp); 

Y(2:NO,2:NO)  =  rot90(temp,2); 

Y(NO:N,2:NO)  =  fliplr(temp); 

Y=Y./(max(max(Y)));  %Normalize  the  Gaussian  distribution  to  one 
%%%  End  of  program  %%% 


CIRCULARLY  TRUNCATED  BESSEL  INPUT  EXCITATION 

function  Y  =  crcbess(a,d,N) 

%crcbess.m:  Y  =  crcbess(a,d,N) 

%  Program  for  generating  circular  Bessel  excitation  functions. 

%a  is  the  width  scaling  factor. 
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%d  is  the  diameter  of  the  circle.  (ODD  inteq~r) 

%N  is  the  width  of  the  square  base.  (EVEN  integer) 

%Example:  z  =  crcbess(l,33,64); 

%  Adopted  from  [Ref  3] 

%Check  that  d  is  an  odd  integer, 
if  rem(d,2)  <0.1; 

error('The  diameter  of  the  circle  must  be  an  ODD  integer'); 

else; 

end; 

%Check  that  N  is  an  even  integer, 
if  rem(N,2)  0.0; 

error('The  width  of  the  square  base  must  be  an  EVEN  integer'); 

else; 

end; 

NO  =  (N/2)+l ;  %NO  is  the  center  location 
r  =  d/2;  %r  is  the  radius  of  the  circle 

Y  =  zeros(N); 
temp  =  zeros(NO-l); 
for  m  =  l:r+l; 

for  n=  l:r+l; 

X  =  sqrt((m-l)^2  +  (n-l)^2); 
if  X  <=  r; 

temp(m,n)=  bessel(0,a*x); 

end; 

end; 

end; 

Y(NO:N,NO:N)  =  temp; 

Y(2:NO,NO:N)  =  flipud(temp); 
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Y(2;NO,2:NO)  =  rot90(temp,2); 
Y(N0:N,2:N0)  =  fliplr(temp); 
%%%  End  of  program  %%% 


APPENDIX  B.  SOURCE  CODE  FOR  FILTER  FUNCTION 
The  followings  are  the  MATLAB  source  code  for  computing  the  propagation 
spatial  filter  function. 

PROPAGATION  SPATIAL  FILTER  FUNCTION 


%%  ioptfil.m 

%%  This  program  generates  the  Bessel  filter  function 
%%  related  files/variables  :  optvar.mat,  pJ164x,  PROPl 
%  Written  by  Nicholas  Lee,  Jul  1998 
clear  all; 


!del  pl28Jl*x*.mat 
!del  optvar28.mat 
N  =  128; 

M=  128; 

NO  =  (N/2)+l; 

Step  =  3; 
z=  lOOe-3; 
time_max  =  .95e-9; 
rho  =  200; 
c  =  3e8; 


%  delete  old  data  files 
%  delete  old  data  files 
%  size  of  square  array 
%  number  of  time  slices 
%  defines  center  of  the  square  array 
%  number  of  leading  zero  time  slices 
%  distance  to  the  observation  plane 
%  time  at  the  final  time  slice 
%  spatial  radius  of  the  filter[sqrt(rhox^2+  rhoy'^2)] 
%  velocity  of  the  light  wave 


%%  Initialize  matrices  to  save  processing  time 
PROP  =  zeros(NO); 

PROPl  =  zeros(N); 

temp  =  zeros(NO);  %bessel  function  of  order  one,  J1 

arg  =  zeros(NO); 

rho_m  =  zeros(NO,l); 

row  =  zeros(NO); 

time  =  zeros(M-Step,l); 
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%%  Generate  M-Step  time  slices  between  z/c  and  time_max. 
time  =  linspace(z/c,  time_max,  M-Step); 

%%  Generate  NO-1  values  of  rho_m  from  0  to  rho. 
rhom  =  linspace(0,rho,NO-l); 

%%  Add  additional  increment  to  rho_m  to  compensate  for  the  off-center 
%%  orientation  of  the  final  NXN  matrix 

rho_m  =  [rho_m  (rho_m(NO-l)+rho_m(2))];  %use  2  b'cos  Matlab  indexing 

%  start  at  l,2,etc 

%  Create  two  NO  x  NO  arrays  of  rho  values  for  function  evaluation. 

[rhox,rhoy]  =  meshgrid(rho_m,rho_m); 

%%  Calculate  matrix  of  radial  distance  values  outside  the  loop 
row  =  sqrt(rhox.^2  +  rhoy.^2); 

%%  Save  those  variables  necessary  for  ioptprop.m  in  a  data  file  optvar.mat 
save  optvar  N  M  NO  Step  time  c  z  row; 

MM=movie(125); 

%%%START  LOOP%%% 
for  m  =  1  :M-Step 

fprintf(  '%3.0f  ,m);  %show  m  value  on  screen 
%Generate  PROP  matrices  Avith  suffix  of  "A"  corresponding 
%  to  the  values  of  the  time  vector 

%  Create  an  NO  x  NO  array  of  argument  values  for  the  bessel  function 
if  time(m)=z/c  %creat  t=z/c  term 

PROP=flipud(2/c-(z.*row.^2)); 

PROPl(l  ;NO,l  :NO)  =  fliplr(PROP); 

PROPl(l:NO,NO:N)  =  PROP(l:NO,l:NO-l); 

PROPl(NO:N,l  :N)  =  flipud(PROPl(2:NO,l  :N)); 

else 

sq  =  sqrt(  c^2*(time(m))^2-z^2  ); 
arg  =  row*sq; 
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%  Evaluate  row*J_l  at  each  argument  value;  create  an  NO  x  NO  array 
temp  =  flipud((-2*z)*(row/sq).*besselj(l,arg)); 

PROP  1(1  :NO,  1  :NO)  =  fliplr(temp); 

PR0P1(1  :NO,NO:N)  =  temp(l  :NO,l  :NO-l); 

PROPl(NO:N,l:N)  =  flipud(PROPl(2:NO,l:N)); 
end 

mesh(PROPl) 

grid  on;  xlabel('fx');  ylabel('fy');  zlabel('amplitude'); 

MM(:,m)=getframe; 

%Correlate  the  name  of  the  variable  PROP  with  the  time  index;ie,  PROPl,  PROP2  etc 
vname  =  ['PROP  1  ',int2str(m)] ;  %set  up  name 

eval([vname,'=  PROPl  ;']); 

%Save  applicable  PROP  in  a  file  named  pJl(N)x(m)A;.e.g.,  PROP15A  in  pJ164x5A 
eval(['save  pJl',int2str(N),'x',int2str(m),' vname]); 
eval(['clear  PROPl  vname]); 
end 

%%%END  LOOP%%% 
movie  (MM); 

%%%  End  of  program  %%% 
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APPENDIX  C.  SOURCE  CODE  FOR  OUTPUT  FIELD  COMPUTATION 
The  followings  are  the  MATLAB  source  code  for  computing  the  output  field 
distribution. 

OUTPUT  FIELD  COMPUTATION 

%  ioptprop.m 

%performs  transient  optical  wave  propagation  simulations 
%It  uses  the  NXN  arrays  "p(N)x(m)A/B"  to 
%  compute  the  propagation  transfer  function. 

%  Size  of  the  variables  NXN  -  input  functions;  M-Step  -  time  slices. 

%  NxM-output_plot 

%  circular,  square  and  gaussian  excitation 
%  Written  by  Nicholas  Lee,  Jul  1998 
clear  all; 

!  del  opt*  .met  %  delete  old  data  files 

%  Load  the  defining  parameters  specified  in  OPTFIL.m 
load  optvar28.mat 

%  Generate  the  INPUT  function;  plot  it. 

N 

disp('N  is  the  width  of  the  base  for  each  function') 
dispC ') 

disp('Please  select  the  excitation  function') 

dispC  1  -  Circle  ') 

dispC  2  -  Table  ') 

dispC  3  -  Circular  Gaussian  ') 

dispC  4  -  Circular  Bessel  ') 

dispC  ’) 

dispC  Please  strike  "Enter"  after  selection.') 
dispC  ’) 

dispC  Default  values  are  in  [  ].') 
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input_func  =  input('Please  enter  an  excitation  function  number  [1] '); 
if  isempty(input_func) 

input_fimc  =  1 
end 

if  input_func  =  1 , 

d  =  input('Please  enter  ODD  diameter,  [25],  d  =  '); 
if  isempty(d) 
d  =  25 
end 

shft_input  =  crcle(d,N); 
elseif  input_func  ==  2, 

w  =  input('Please  enter  ODD  width,  [25],  w  = '); 
if  isempty(w) 
w  =  25 
end 

shft_input  =  table(w,N); 

d  =  w; 

elseif  input_func  =  3, 

sigma  =  input('Please  enter  the  standard  deviation,  [12], sigma  =  '); 
if  isempty(sigma) 
sigma  =  12 
end 

d  =  input('Please  enter  the  ODD  diameter,  [25],  d  = '); 
if  isempty(d) 
d  =  25 
end 

shft_input  =  crcgaus(sigma,d,N); 
elseif  inputfunc  ==  4, 

a  =  input('Please  enter  the  width  scaling  factor,,  [1],  a= '); 
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if  isempty(a) 
a=l 
end 

d  =  input('Please  enter  the  ODD  diameter,  [25],  d  = '); 
if  isempty(d) 
d  =  25 
end 

shfl_input  =  crcbess(a,d,N); 
else 

error('Incorrect  Excitation  Function  Selection') 
end 

%%  Shift  input  quadrants  and  take  the  2-D  FFT  to  produce  F_INPUT. 
input  =  (fftshift(shft_input)); 

F_input  =  real(fft2(input)); 

%  Shift  F_input  in  preparation  of  multiplication  with  PROPl 
Fshft_input  =  fftshift(F_input); 

%  Array-multiply  the  filter  transfer  function  PROPl  and  Fshft  input. 
dispCPerforming  array  multiplication....'); 

%%%  Start  loop  %%% 
for  m  =  1  :M-Step 

Q)rintf(  '%2.0f,  ',m) 
pause(l) 

%%  Load  filter  transfer  function 

filenamel  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 

%  Array-multiply  filter  transfer  function  with  Fshft_input 
Fshftoutputl  =  vnamel.*(Fshft_input); 

%Clear  unnecessary  variables  to  free  RAM 
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clear  vnamel; 

eval(['clear  PROP  1  ',int2str(in), 

%  Shift  Fshft_outputl  to  comer  geometry  prior  to  taking  the  IFFT2 
F_outputl  =  fftshift(Fshft_outputl); 

%  Take  IFFT  of  F_outputl 

output=ifft2(F_output  1 ) ; 

%  Shift  output  1  prior  to  summation««««OUTPUT 
shft_output  =  fftshift(output); 

%View  the  magnitude  of  the  shifted  output<««INTENSITY 
shft_outabs  =  abs(shft_output); 
shft_intensity  =  (shft_outabs).^2; 

%Shft_outabs  as  outabs  and  store  into  file  optab(time(m)) 
vname  =  [' outabs', int2str(m)]; 
eval([vname,'=shft_outabs ;']) 
eval(['save  optab',int2str(m),'  ',vname]) 
vname  =  ['inten',int2str(m)]; 
eval([vname,'=shft_intensity ;']) 
eval(['save  optint',int2str(m),' vname]) 

%Save  the  NO  row  of  the  magnitude  of  the  shifted  output  in  the 
%m+Step  column  of  output_plot. 

output_plot(  1  :N,m+Step)=shft_outabs(NO,  1  :N)'; 
end 

%%%%  End  loop  %%% 

%Save  contents  of  "output_plot"  as  NxM  array, 
filename  =  ['op',int2str(d),'x',int2str(M)];  %  File:  op(d)x(M) 
eval(['save  ',filename,'  output_plot']); 

%  plot  the  responses  at  each  stage 

save  excit  shftinput  input  F_input  Fshftinput  output jplot 
%%%  End  of  program  %%% 
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APPENDIX  D.  SOURCE  CODE  FOR  2D  AND  3D  GRAPHICS 


The  followings  are  the  MATLAB  source  code  for  plotting  all  two-  and  three- 
dimensional  graphics. 

PLOT  FILTER  FUNCTION 
%plotfilter.m 

%This  program  plots  all  the  filter  function 
%It  uses  data  files  pj  1 128xn.mat 

%Written  by  Nicholas  Lee,  Jul  1998 

clear  all 

cs60  =  zeros(128); 

csl25  =  zeros(128); 

load  pj  11 28x1. mat 
load  pj  1 128x2.mat 
load  pj  1 128x4.mat 
loadpjll28x8.mat 
load  pj  1 128x60.mat 
load  pjl  128xl25.mat 
figure(l) 
subplot(  1,2,1) 
mesh(PROPl  1) 

axis  square,  title('(a)  Time  slice  1') 

grid  on,  xlabel(’fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -9000  0]) 

set(gca,'xtick’,[0,65,128],'ytick’,[0,65,128],'ztick',[-9000,-6000, -3000,0]) 
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subplot(  1,2,2) 
mesh(PR0P12) 

axis  square,  title('(b)  Time  slice  2') 

grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -900  1200]) 

set(gca,'xtick',[0,65, 1 28]  ,'ytick',  [0,65, 1 28], 'ztick',  [-800,-400,0,400,800, 1 200]) 

figure(2) 
subplot(  1,2,1) 
mesh(PROP14) 

axis  square,  title('(d)  Time  slice  4') 

grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -600  600]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],’ztick’,[-600,-400,-200,0,200,400,600]) 

subplot(l  ,2,2) 
mesh(PROP18) 

axis  square,  title('(d)  Time  slice  8') 

grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -600  600]) 

set(gca,'xtick',[0,65,128],'ytick’,[0,65,128],’ztick’,[-600,-400,-200,0,200, 400,600]) 

figure(3) 
subplot(  1,2,1) 
mesh(PROP160) 

axis  square,title('(a)  Time  slice  60') 

grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -50  50]) 
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set(gca,'xtick', [0,65, 128], 'ytick',[0, 65, 128], 'ztick', [-40, -20, 0,20,40]) 


subplot(l,2,2) 

mesh(PROP1125) 

axis  square,  title('(b)  Time  slice  125') 

grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -20  20]) 

set(gca,'xtick',  [0,65,128]  ,'ytick',  [0,65,1 28]  ,'ztick',  [-20,- 1 0,0, 1 0,20]) 

cs60(65,l:128)  =  PROP160(65,l:128); 

csl25(65,l:128)=PROPl  125(65,1:128); 

figure(4) 

subplot(  1,2,1) 

mesh(cs60) 

axis  square, title('(c)  Time  slice  60') 

grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -50  50]) 

set(gca,'xtick',[0,65,128],'ytick',[0, 65,128], 'ztick',[-40,-20,0, 20,40]) 
view(0,0) 

subplot(  1,2,2) 
mesh(csl25) 

axis  square,  title('(d)  Time  slice  125') 

grid  on,  xlabel(’fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

axis([0  128  0  128  -20  20]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[-20,-10,0,10,20]) 

view(0,0) 

%%%  End  of  program  %%% 
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PLOT  FIELD 


%plotfield.m 

%This  program  plots  all  the  graphics  for  the  input  and  output  fields 
%It  uses  data  files  excit,  optvar.mat,  optabm.mat  and  optintm.mat 

%  Written  by  Nicholas  Lee,  Jul  1998 
clear  all 

%  Load  the  defining  parameters  specified  in  OPTFIL.m 
load  excit.mat 

%load  output  field 

load  optabl.mat 

load  optabSO.mat 

load  optablOO.mat 

load  optabl25.mat 

%load  intensity  to  create  image 

load  optintl  .mat 

load  optint50.mat 

load  optintl  OO.mat 

load  optintl25.mat 

figure(l)  %  input  excitation 

subplot(  1,2,1) 

mesh(shft_input);titleC(a)SHFT-INPUT'); 
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axis  square; 
axis([0  128  0  128  0  1]) 

set(gca,'xtick',[0,65,128];ytick',[0,65,128];ztick',[0,0.5,l]) 
grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Z-axis') 

subplot(  1,2,2) 

mesh(input);title(’(b)  INPUT') 

axis  square; 

axis([0  128  0  128  0  1]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick', [0,0.5,!]) 
grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Z-axis') 
pause;  ' 

figure(2)  %  Fourier  transform  of  shifted  input 
subplot(  1,2,1) 

mesh(F_input);title('(c)  F-INPUT ') 
axis  square; 

axis([0  128  0  128  -100  500]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[-100,0,100,200,300,400,500]) 
grid  on,  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

subplot(  1,2,2) 

mesh(Fshft_input);title('(d)  FSHFT-INPUT ') 
axis  square; 

axis([0  128  0  128  -100  500]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[-100,0,100,200,300,400,500]) 

grid  on,  xlabel('f!C-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 

pause; 
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figvire(3) 

subplot(2,2,l) 

mesh(outabsl) 

axis  square, title('(a)  Output  field  at  time  slice  1') 

grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel(' Amplitude') 

axis([0  128  0  128  0  1500]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],’ztick’,[0,500,1000,1500]) 

subplot(2,2,2) 

mesh(intenl) 

axis  square,title('Image  at  time  slice  1') 

grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Z-axis') 

axis([0  128  0  128  0  2e6]) 

colorbar 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[  ]) 
view(0,90) 

subplot(2,2,3) 

mesh(outabs50) 

axis  square,title('(b)  Output  field  at  time  slice  50') 

grid  on,  xlabel('X-axis'),  ylabelf' Y-axis'),  zlabelf' Amplitude') 

axis([0  128  0  128  0  30]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick’,[0,10,20,30]) 

subplot(2,2,4) 

mesh(inten50) 

axis  square,title('Image  at  time  slice  50') 

grid  on,  xlabel('X-axis'),  ylabel('Y-axis'),  zlabel('Z-axis') 

axis([0  128  0  128  0  600]) 
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colorbar 

set(gca,'xtick', [0,65, 128], 'ytick', [0,65, 128], 'ztick',[  ]) 

view(0,90) 

pause; 

figure(4) 

subplot(2,2,l) 

mesh(outabslOO) 

axis  square, title('(c)  Output  field  at  time  slice  100') 

grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel(' Amplitude') 

axis([0  128  0  128  0  20]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65, 128], 'ztick',[0,5, 10,15,20]) 

subplot(2,2,2) 

mesh(intenlOO) 

axis  square,title('Image  at  time  slice  100') 

grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Z-axis') 

axis([0  128  0  128  0  250]) 

colorbar 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[  ]) 
view(0,90) 

subplot(2,2,3) 

mesh(outabsl25) 

axis  square,title('(d)  Output  field  at  time  slice  125') 

grid  on,  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel(' Amplitude') 

axis([0  128  0  128  0  10]) 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[0,2,4,6,8,10]) 
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subplot(2,2,4) 

mesh(intenl25) 

axis  square, title('Image  at  time  slice  125') 

grid  on,  xlabel('X-axis'),  ylabel('Y-axis'),  zlabel('Z-axis') 

axis([0  128  0  128  0  80]) 

colorbar 

set(gca,'xtick',[0,65,128],'ytick',[0,65,128],'ztick',[  ]) 

view(0,90) 

pause; 

figure(5)  %  output_plot  with  close-up  view 
subplot(l,2,l) 

mesh(output_plot);title('(a)  Total  output'); 
axis  square; 

axis([0  128  0  128  0  800]) 

set(gca,'xtick',[0,40,80,128],’ytick’,[0,65,128],'ztick',[0,200,400,600,800]) 
grid  on,  xlabel('time  slice'),  ylabel('space'),  zlabel('amplitude') 
view(52.5,30) 

subplot(  1,2,2) 

mesh(output_plot);title('(b)  Magnified  view  of  total  output'); 
axis  square; 

axis([0  128  0  128  0  80]) 

set(gca,’xtick',[0,40,80,128],’ytick',[0,65,128],’ztick’,[0,20,40,60,80]) 
grid  on,  xlabel('time  slice'),  ylabel('space'),  zlabel('amplitude') 
view(52.5,30) 
pause; 

figure(6)  %  output_plot  side  &  front  views 
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subplot(  1,2,1) 

mesh(output_plot);title('(a)  Close-up  cross-section  view  of  total  output'); 
axis  square; 

axis([0  128  0  128  0  80]) 

set(gca,'xtick',[0,20, 40,60, 80, 100,128],'ytick',[0,65,128],'ztick',[0, 20,40,60,80]) 

grid  on,  xlabel('time  slice'),  ylabel('space'),  zlabel('amplitude ') 

view(0,0) 

subplot(  1,2,2) 

mesh(ou1put_plot);title('(b)  Close-up  front  view  of  total  output'); 
axis  square; 

axis([0  1-28  0  128  0  80]) 

set(gca,'xtick',[0,40,80,128],'ytick',[0,30,65,90,128],'ztick',[]) 
grid  on,  xlabel('time  slice'),  ylabel('space'),  zlabel(' ') 
view(90,45) 

%%%  End  of  program  %%% 
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APPENDIX  E.  SOURCE  CODE  FOR  ANIMATION  PROGRAMS 
The  followings  are  the  MATLAB  source  code  for  the  animation  programs. 


ANIMATION  FORMAT  1 
%  animate l.m 

%This  program  animate  filter  function ,  output  field  and  image 
%uses  N=64 

%Written  by  Nicholas  Lee,  Aug  1998 
clear  all; 

%  Load  the  defining  parameters  specified  in  IOPTFIL.m 
load  optvar.mat 

%  Array-multiply  the  shifted  transfer  function  PRROP  and  Fshft_input. 
dispC Animation  in-progress....'); 

movie_figure  =  figure('position',[50  200  600  220]);%col,  row 
MM=moviein(M-Step,movie_figure); 

%%%  Start  loop  %%% 
for  m  =  1  :M-Step 

nic=['time-slice  ',int2str(m)  ]; 

Q)rintf(  '%2.0f,  ',m); 
ifm  =1 

filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval([’vnamel=PROPr,int2str(m), 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m), 
subplot(  1,3,1) 

mesh(vnamel);title('(a)  Filter  spatial  fi'equency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
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axis([0  64  0  64  -9000  0]) 

set(gca,'xtick',  [0,3 3 ,64],'ytick',  [0,3 3 ,64],'ztick',  [-9000,0]) 
subplot(l,3,2) 

mesh(vname2);title('  (b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  1500]) 

set(gca,'xtick', [0,33,64], 'ytick',[0,33,64],'ztick’,[0,500, 1000, 1500]) 
subplot(l,3,3) 

mesh(vname2);title('(c)  Image’) 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  1500]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0, 33,64], 'ztick',[0,500, 1000,1500]) 
view(0,90) 

eval(['text(15,-30,0,nic);']); 
elseif  m  <=  3 

filename  1  =  ['pJr,int2str(N),'x’,int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ’,int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(l,3,l) 

mesh(vnamel);title('(a)  Filter  spati^ll  frequency  response') 
axis  square; 
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grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -400  600]) 

set(gca;xtick’,[0, 33, 64];ytick', [0,33, 64], 'ztick', [-400, -200, 0,200, 400,600]) 
subplot(l,3,2) 

mesh(vname2);title('  (b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  250]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,50,l  00, 1 50,200,250]) 
subplot(l,3,3) 

mesh(vname2);title('(c)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  250]) 

colorbar; 

set(gca,'xtick’,[0,33,64],'ytick',[0,33,64],'ztick',[0,50, 1 00, 1 50,200,250]) 
view(0,90) 

eval(['text(l  5,-30,0,nic);']); 
elseif  m  <=1 1 

filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(  1,3,1) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
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axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel(Ty-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -300  300]) 

set(gca,'xtick’,[0,33,64],'ytick',[0,33,64],'ztick',[-300,-200,- 

100,0,100,200,300]) 

subplot(l,3,2) 

mesh(vname2);title('  (b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],’ztick',[0,20,40,60,80,100]) 

subplot(l,3,3) 

mesh(vname2);title('(c)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  100]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,20,40, 60,80,100]) 
view(0,90) 

eval(['text(l  5,-30,0,nic);']); 
elseif  m  <=25 

filenamel  =  ['pJr,int2str(N),'x',int2str(m)  ]; 

eval(['load  ',filenamel]); 

eval(['vnamel  =PROP  1  ',int2str(m),';']); 

filename2  =  ['optab',int2str(m)  ]; 

eval(['load  ',filename2]); 

eval(  [' vname2=outabs'  ,int2str(m),';'] ) ; 
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subplot(  1,3,1) 

mesh(vnainel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('j&c-axis'),  ylabel('ly-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -100  100]) 

set(gca,'xtick',  [0,33 ,64]  ,'ytick',  [0,3 3 ,64]  ,'ztick',  [- 1 00,-5 0,0,50, 1 00]) 
subplot(l,3,2) 

mesh(vname2);title('  (b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  50]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,10,20,30,40,50]) 

subplot(l,3,3) 

mesh(vname2);title('(c)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  50]) 

colorbar 

set(gca,'xtick',  [0,3 3 ,64]  ,'ytick',  [0,33 ,64]  ,'ztick',  [0, 1 0,20,3 0,40,50]) 
view(0,90) 

eval(['text(l  5,-30,0,nic);']); 
else  m  <=61 

filename  1  =  ['pJl',int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
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eval(['vname2=outabs',int2str(m), 

subplot(l,3,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -50  50]) 

set(gca,'xtick', [0,33, 64], 'ytick', [0,33, 64], 'ztick',[-50, 0,50]) 
subplot(l,3,2) 

mesh(vnanie2);title('  (b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  30]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick', [0,10,20,30]) 
subplot(l,3,3) 

mesh(vname2);title('(c)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  30]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0, 10,20,30]) 
view(0,90) 

eval(['text(15,-30,0,nic);']); 

end 

figure(mo  vie_figure) ; 

MM( :  ,m)=getframe(gcf) ; 
end 

%%%  END  OF  LOOP  %%% 
echo  off 
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dispC '); 

disp('Press  a  key  to  play  back  movie.'); 

pause 

echo  on 

start_fi-ame=input('Enter  start  frame:'); 
end_frame=mput('Enter  end  frame:'); 
movie(movie_figure,MM,[l  (start_frame:end_frame)],  1); 
echo  off 

%%%End  of  program  %%% 


ANIMATION  FORMAT  2 
%  animate2.m 

%This  program  animate  filter  function ,  output  field,  total  output  (side  view) 
%and  image 
%uses  N=64 

%Written  by  Nicholas  Lee,  Aug  1998 
clear  all; 

%  Load  the  defining  parameters  specified  in  IOPTFIL.m 

load  optvar.mat 

center=zeros(N); 

%  Array-multiply  the  shifted  transfer  function  PRROP  and  Fshft_input. 
dispC Animation  in-progress....'); 

movie_figure  =  figure('position',[50  100  450  350]);%col,  row 
MM=moviein(M-Step,movie_figure); 

%%%  Start  loop  %%% 
for  m  =  1,:M-Step 

nic=['time-slice  ',int2str(m)  ]; 
fprintf(  '%2.0f,  ',m); 
ifm==l 
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filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 

eval(['load ',filenamel]); 

eval(['vname  1  =PROP  1  ',int2str(m), 

filename2  =  ['optab',int2str(m)  ]; 

eval(['load  ',filename2]); 

eval(  [' vname2=outabs'  ,mt2str(m),' ; 

subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -9000  0]) 

set(gca,'xtick', [0,33, 64], 'ytick',[0,33, 64], 'ztick', [-9000,0]) 
subplot(2,2,2) 

mesh(vname2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  1500]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],’ztick',[0,500, 1 000,1 500]) 
subplot(2,2,3) 

center(m,  1  :N)=vname2(NO,  1  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,20,40,60,80,100]) 

view(90,0) 

hold  on 
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subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  1500]) 

colorbar 

set(gca,'xtick', [0,33, 64],'ytick',[0,33, 64], 'ztick', [0,500, 1000, 1500]) 
view(0,90) 

eval(['text(64,- 1 3,0,nic);']); 
elseif  m  <=  3 

filenamel  =  ['pjr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -400  600]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[-400,-200,0,200,400,600]) 

subplot(2,2,2) 

mesh(vn2ime2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  250]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,50, 100,150,200, 250]) 
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subplot(2,2,3) 

center(m,  1  ;N)=vname2(N0, 1  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca;xtick',[0, 33,64], 'ytick’,[0, 33,64], 'ztick', [0,20, 40, 60, 80, 100]) 

view(90,0) 

hold  on 

subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel ('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  250]) 

colorbar; 

set(gca,'xtick’,[0,33,64],'ytick',[0,33,64],'ztick',[0,50,100,150,200,250]) 

view(0,90) 

eval(['text(64,-13,0,nic);']); 
elseif  m  <=1 1 

filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),’;']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  fi'equency  response') 
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axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -300  300]) 

set(gca;xtick',[0,33,64];ytick',[0,33,64],'ztick’,[-300,-200,- 

100,0,100,200,300]) 

subplot(2,2,2) 

mesh(vname2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,20,40,60,80,100]) 

subplot(2,2,3) 

center(m,  1  :N)=vname2(NO,  1  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,3 3 ,64],'ytick',[0,3 3 ,64],'ztick',[0,20, 40,60,80, 1 00]) 

view(90,0) 

hold  on 

subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  100]) 

colorbar 
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set(gca, 'xtick', [0,33, 64], 'ytick', [0,33, 64], 'ztick', [0,20, 40, 60,80, 100]) 
view(0,90) 

eval(['text(64,-13,0,nic);']); 
elseif  m  <=25 

filename  1  =  ['pJr,int2str(N),'x’,int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vnamel  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64-100  100]) 

set(gca,’xtick’,[0,33,64],'ytick',[0,33,64],'ztick',[-100,-50, 0,50,100]) 
subplot(2,2,2) 

mesh(vname2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  50]) 

set(gca,’xtick’,[0,33,64],'ytick',[0,33,64],'ztick',[0,10,20,30,40,50]) 

subplot(2,2,3) 

center(m,  1  :N)=vname2(NO,l  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
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axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64];ytick',[0,33,64],’ztick',[0,20,40,60,80,100]) 

view(90,0) 

hold  on 

subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  50]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,10,20,30,40,50]) 

view(0,90) 

eval(['text(64,- 1 3 ,0,nic);']) ; 
else  m  <=61 

filenamel  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vnamel=PROP  1  ',int2str(m), ';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m), ';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -50  50]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[-50,0,50]) 

subplot(2,2,2) 
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mesh(vnaine2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel('Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  30]) 

set(gca,'xtick', [0,33, 64], 'ytick',[0,33,64],'ztick’, [0,10,20, 30]) 
subplot(2,2,3) 

center(m,  1  :N)=vname2(NO,  1  :N); 
mesb(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0, 20,40,60,80,100]) 

view(90,0) 

bold  on 

subplot(2,2,4) 

mesb(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  30]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0, 10,20,30]) 
view(0,90) 

eval(['text(64,-13,0,nic);']); 

end 

figure(movie_figure); 
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MM(:,m)=getfranie(gcf); 


end 

%%%  END  OF  LOOP  %%% 
echo  off 

dispC); 

dispCPress  a  key  to  play  back  movie.'); 

pause 

echo  on 

start_frame=input('Enter  start  frame:'); 
end _frame=input('Enter  end  frame:'); 
movie(movie_figure,MM,[l  (start_frame:end_frame)],l); 
echo  off- 

%%%  End  of  program  %%% 


ANIMATION  FORMAT  3 

%  animates  .m 

%This  program  animate  filter  function ,  output  field,  total  output  (magnified) 
%and  image 
%uses  N=64 

%Written  by  Nicholas  Lee,  Aug  1998 
clear  all; 

%  Load  the  defining  parameters  specified  in  lOPTFIL.m 

load  optvar.mat 

center=zeros(N); 

%  Array-multiply  the  shifted  transfer  function  PRROP  and  Fshft_input. 
dispC Animation  in-progress....'); 

movie_figure  =  figure('position',[50  100  450  350]);%col,  row 
MM=moviein(M-Step,movie_figure); 

%%%  Start  loop  %%% 
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for  m  =  1  :M-Step 

nic=['time-slice  ',int2str(m)  ]; 
fprintf(  '%2.0f,  ',m); 
ifm=l 

filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
evaI(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m), 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -9000  0]) 

set(gca,'xtick',[0, 33,64], ’ytick',[0, 33,64], 'ztick', [-9000,0]) 
subplot(2,2,2) 

mesh(vname2);title(’(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  1500]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,500,1000,1500]) 

subplot(2,2,3) 

center(m,l  :N)=vname2(NO,l  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
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axis([0  64  0  64  0  100]) 

set(gca,'xtick', [0,33, 64], 'ytick',[0, 33,64], 'ztick',[0, 20,40, 60, 80, 100]) 

view(142.5,30) 

hold  on 

subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  1500]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,500,1000,1500]) 

view(0,90) 

eval(['text(64,-13,0,nic);']); 
elseif  m  <=  3 

filename  1  =  ['pjr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
ev£il(['vname2=outabs',int2str(m),';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -400  600]) 

set(gca,'xtick’,[0,33,64],'ytick',[0,33,64],'ztick',[-400,-200,0, 200,400,600]) 
subplot(2,2,2) 
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mesh(vname2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  250]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,50,100,150,200,250]) 

subplot(2,2,3) 

center(m,  1  :N)=vname2(NO,  1  ;N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca, 'xtick', [0,33, 64], 'ytick', [0,33, 64], 'ztick', [0,20, 40, 60, 80,100]) 
view(  142.5,30) 
hold  on 

subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  250]) 

colorbar; 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,50,100,150,200,250]) 

view(0,90) 

eval(['text(64,-13,0,nic);']); 
elseif  m  <=1 1 

filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vnamel  =PROP  1  ',int2str(m),';']); 
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filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m), 

subplot(2,2,l) 

mesh(vnainel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fic-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64  -300  300]) 

set(gca;xtick',[0,33,64],'ytick',[0,33,64],’ztick',[-300,-200,- 

100,0,100,200,300]) 

subplot(2,2,2) 

mesh(vname2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33 ,64],'ytick',[0,33,64],'ztick',[0,20,40,60,80, 1 00]) 
subplot(2,2,3) 

center(m,  1  ;N)=vname2(NO,  1  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,20,40,60,80,100]) 

view(142.5,30) 

hold  on 

subplot(2,2,4) 
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mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis’),  ylabel('Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  100]) 

colorbar 

set(gca,'xtick',[0, 33, 64], 'ytick’,[0, 33,64], 'ztick', [0,20, 40, 60,80, 100]) 
view(0,90) 

eval(['text(64,-l  3,0,nic);']); 
elseif  m  <=25 

filename  1  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load ',filenamel]); 
eval(['vname  1  =PROP  1  ',int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
axis([0  64  0  64-100100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[- 100,-50,0,50, 100]) 
subplot(2,2,2) 

mesh(vname2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 
axis([0  64  0  64  0  50]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,10,20, 30,40,50]) 
subplot(2,2,3) 
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center(m,l  ;N)=vnajtne2(N0, 1  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Tinie-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick’,[0,33,64],'ytick’,[0,33,64],'ztick',[0,20,40,60,80,100]) 

view(142.5,30) 

hold  on 

subplot(2,2,4) 

mesh(vnanie2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  50]) 

colorbar 

set(gca,'xtick', [0,33, 64], 'ytick', [0,33, 64],'ztick', [0,10, 20,30, 40, 50]) 
view(0,90) 

eval(['text(64,-13,0,nic);']); 
else  m  <=61 

filenamel  =  ['pJr,int2str(N),'x',int2str(m)  ]; 
eval(['load  ',filenamel]); 
eval(['vnamel=PROPr,int2str(m),';']); 
filename2  =  ['optab',int2str(m)  ]; 
eval(['load  ',filename2]); 
eval(['vname2=outabs',int2str(m),';']); 
subplot(2,2,l) 

mesh(vnamel);title('(a)  Filter  spatial  frequency  response') 
axis  square; 

grid  on;  xlabel('fx-axis'),  ylabel('fy-axis'),  zlabel('fz-axis') 
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axis([0  64  0  64  -50  50]) 

set(gca,'xtick',[0,33,64],'ytick’,[0,33,64],'ztick',[-50,0,50]) 

subplot(2,2,2) 

mesh(vnaine2);title('(b)  Image  field  intensity') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel ('Intensity') 
axis([0  64  0  64  0  30]) 

set(gca, 'xtick', [0,33, 64], 'ytick',[0, 33,64], 'ztick', [0,10, 20, 30]) 
subpIot(2,2,3) 

center(m,  1  ;N)=vname2(N  0, 1  :N); 
mesh(center);title('(c)  Field  distribution') 
axis  square; 

grid  on;  xlabel('Space'),  ylabel('Time-slice'),  zlabel('Intensity') 
axis([0  64  0  64  0  100]) 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,20,40, 60,80,100]) 
vie  w( 142. 5, 30) 
hold  on 

subplot(2,2,4) 

mesh(vname2);title('(d)  Image') 
axis  square; 

grid  on;  xlabel('X-axis'),  ylabel(' Y-axis'),  zlabel('Intensity') 

axis([0  64  0  64  0  30]) 

colorbar 

set(gca,'xtick',[0,33,64],'ytick',[0,33,64],'ztick',[0,10,20,30]) 

view(0,90) 

eval(['text(64,-13,0,nic);']); 


no 


end 

figure(movie_figiire); 

MM(:,m)=getframe(gcf); 

end 

%%%  END  LOOP  %%% 
echo  off 
dispC '); 

disp('Press  a  key  to  play  back  movie.'); 

pause 

echo  on 

start_fi:ame=input('Enter  start  frame:'); 
end_frame=input('Enter  end  frame:'); 
movie(movie_figure,MM,  [  1  (start_frame:end_frame)],  1 ); 
echo  off 

%%%  End  of  program  %%% 
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